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1  Introduction 


During  this  task,  we  established  a  method  for  extracting  the  most  significant  basis  functions 
from  a  physics-based  model  of  the  thermosphere  and  ionosphere.  A  new  model  was  then  con¬ 
structed  by  replacing  the  spherical  harmonic  functions  used  in  an  existing  semi-empirical  ther¬ 
mosphere  model  with  these  basis  functions.  The  resulting  hybrid  model,  published  in  the  Space 
Weather  Journal  ( Sutton  etal.  |[7lJ|),  demonstrated  the  ability  to  represent  sparse  thermospheric 
data  in  a  more  accurate  and  efficient  way.  In  addition,  we  chose  Jacchia  1970  as  our  underlying 
semi-empirical  model  so  that  the  new  hybrid  model  functions  as  a  direct  replacement  for  the  High 
Accuracy  Satellite  Drag  Model  [H|  (HASDM)  used  by  AFSPC  and  the  JSpOC. 

We  also  worked  to  improve  the  hybrid  model  by  investigating  the  limitations  of  and  possible 
improvements  to  the  underlying  assumptions  imposed  by  the  chosen  models.  In  particular,  we 
realized  that  there  was  a  large  and  potentially  significant  component  missing  from  all  general 
circulation  models,  that  being  the  appropriate  treatment  of  helium  in  the  upper  thermosphere. 
While  helium  concentration  in  the  lower  thermosphere  constitutes  merely  a  minor  constituent,  its 
small  atomic  mass  compared  with  that  of  N2,  02,  and  O  causes  its  concentration  to  decrease 
much  more  slowly  with  altitude.  At  some  height,  helium  overtakes  N2,  02,  and  O  as  the  dominant 
species.  Generally  this  occurs  between  400  and  800  km,  depending  on  season,  location,  and 
solar  flux.  At  this  transition  altitude,  the  scale  height — which  gives  an  indication  of  the  exponential 
rate  of  decay  of  mass  density  with  altitude — quickly  transitions  from  ~50  km  to  ~200  km. 

Not  only  is  the  omission  of  helium  an  issue  for  physical  models,  but  the  Jacchia  1970  model, 
JB2006,  and  JB2008  all  suffer  from  this  problem  as  well.  Because  these  models  do  not 
correctly  account  for  helium,  their  vertical  structures  can  be  grossly  inaccurate  in  the  upper 
thermosphere.  The  impact  on  the  operational  system  is  clearly  seen  from  the  well-known  fact 
that  HASDM  quickly  loses  accuracy  at  high  altitudes  and  high  latitudes,  i.e.  where  helium 
concentrations  are  highest.  In  effect,  the  bulk  of  the  data  used  to  drive  HASDM  comes  from 
altitudes  around  300-500  km;  however,  this  information  is  being  inaccurately  extrapolated  to 
higher  altitudes  around  700-1000  km,  a  region  through  which  more  than  10,000  trackable 
objects  travel  (or  about  75%  of  all  trackable  LEO  objects). 

To  rectify  this  situation,  we  undertook  the  project  of  adding  helium  to  the  solution  of  a  popular 
first-principles  model  of  the  thermosphere  and  ionosphere.  In  doing  this,  we  discovered  that  this 
could  only  be  done  self-consistently  by  solving  for  helium  as  a  major  species,  as  opposed  to  as  a 
minor  species.  The  product  of  this  work  was  the  first  physics-based  model  to  accurately  simulate 
the  seasonal  and  latitudinal  behavior  of  the  helium  distribution.  This  work  was  published  in  the 
Journal  of  Geophysical  Research  ( Sutton  eta I.  |JZ2|). 
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2  Thermospheric  Basis  Functions 
2.1  Background 

Atmospheric  drag  is  the  dominant  and  most  difficult  force  to  determine  and  predict,  in  the  orbit 
propagation  model  of  low  earth  orbiting  satellites  [36].  The  drag  acceleration  vector,  a,  acting  on 
a  satellite  can  be  related  to  several  satellite  and  orbital  properties  and  to  the  atmospheric  mass 
density,  p,  by  the  following  equation: 

d=-l/2(CDAref/m)p\V\V  (1) 

where  CD  is  the  coefficient  of  drag,  Aref  is  the  reference  satellite  area  projected  into  the  ram 
direction,  m  is  the  satellite  mass,  and  V  is  the  satellite  velocity  vector  with  respect  to  the  at¬ 
mosphere.  Changes  in  neutral  density  typically  contribute  the  most  to  the  total  variability  of 
drag  acceleration,  however,  the  ( CDAref/m )  term  and  thermospheric  winds  (entering  through 
the  \V\V  term)  can  also  contribute  significant  amounts  at  times. 

The  thermosphere  is  a  strongly  driven  dynamic  system.  Variations  of  neutral  density  in  the 
thermosphere  depend  not  only  on  location  but  on  solar  and  geophysical  conditions  as  well. 
Satellite  accelerations  can  exhibit  changes  of  more  than  an  order  of  magnitude  over  the  course 
of  an  11 -year  solar  cycle,  and  of  a  factor  of  2-4  during  moderate  geomagnetic  events  with  an 
approximate  time-scale  of  1-2  days.  During  geomagnetic  storms,  the  spatial  distribution  and 
temporal  response  strongly  depends  on  latitude  and  local  time,  both  in  geographic  and  magnetic 
coordinates. 

The  Air  Force  Space  Command  (AFSPC)  requires  accurate  knowledge  and  forecasting  of 
neutral  density  between  90  and  500  km.  This  requirement  is  focused  on  improving  the  effi¬ 
ciency  of  satellite  catalog  maintenance,  the  accuracy  of  reentry  predictions,  and  the  reliability  of 
satellite  conjunction  analysis.  In  support  of  this  goal,  the  current  pursuit  focuses  on  the  incre¬ 
mental  improvement  of  calibration  techniques  for  semi-empirical  thermospheric  models.  This  is 
accomplished  by  finding  the  Principal  modes  of  the  thermosphere  as  defined  by  a  physics-based 
model.  This  model-driven  approach  is  taken  to  avoid  the  obstacles  of  insufficient  data  coverage 
that  are  encountered  when  constructing  orthogonal  functions  using  data-driven  techniques  [see 


2.2  Methods,  Assumptions,  and  Procedures 


2.2.1  The  Jacchia  Semi-Empirical  Model 


The  Jacchia  1970  model  (J70  hereafter)  [22]  is  a  static  diffusion  model  of  the  upper  atmo¬ 
sphere.  The  thermospheric  portion  of  the  model  begins  at  1 05  km.  Above  this  level,  J70  specifies 
the  total  number  density,  ni:  of  each  species,  i  =  {Ar,  He,  N2,  O,  or  02},  by  vertically  integrating 
the  diffusion  equation: 


drii 
n  i 


mig 

kBT 


dz  — 


dT 

~T 


(1  +  ati) 


(2) 


where  m  is  the  molecular  weight,  g  is  the  acceleration  due  to  gravity,  2  is  the  height,  kB  is  the 
Boltzmann  constant,  T  is  the  temperature,  and  a  is  the  thermal  diffusion  coefficient  ( a  assumed 
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to  be  -0.38  for  He  and  0  for  all  other  species).  The  total  mass  density,  p,  can  be  calculated  as 
the  summation  of  the  individual  mass  densities  of  each  species:  p  =  Y.inimi/NA,  where  NA  is 
Avogadro’s  number.  The  vertical  temperature  profile,  T(z),  required  to  carry  out  the  integration 
in  (2)  is  parameterized  in  the  following  form: 


(3) 


where  A  =  T, -  Tx),  B  =  4.5xl0-6  km”1,  Tx  is  the  temperature  at  a  prescribed  inflection 
point  zx  =  125  km,  Gx  is  the  temperature  gradient  at  the  inflection  point,  is  the  exospheric 
temperature,  and  j  is  set  to  2.5.  With  the  application  of  several  additional  constraints  [see  22], 
the  model  produces  a  vertical  density  profile  that  is  uniquely  specified  by  Tx  and  T^.  These 
temperatures  have  been  fit  to  satellite  drag  data  and  are  represented  by  parametric  functions  of 
solar  flux,  geomagnetic  activity,  latitude,  and  apparent  solar  local  time. 

Marcos  et  al.  [37]  first  attempted  to  calibrate  this  model  by  estimating  a  global  correction  to 
the  nighttime  minimum  value  of  TE,  designated  as  A Tc,  to  bring  the  model  into  better  agreement 
with  recent  satellite  tracking  data.  Storz  et  al.  [69]  extended  this  technique  by  estimating  a 
spherical  harmonic  (SH)  field  for  A Tc  as  well  as  for  a  correction  to  the  value  of  Tx,  designated  as 
A Tx.  Both  global  temperature  corrections  take  the  following  form: 


oo 


A T(<p,6)  =  EE  PJn( sin  ip)  {Cr  cos  md  +  S sin  m6} 


(4) 


1=0  m= 0 


where  the  P’s  are  orthonormalized  associated  legendre  polynomials  of  degree  and  order  [Ixm], 
the  C’s  and  S’ s  are  the  expansion  coefficients  with  Sf  omitted,  p  is  latitude,  and  9  is  local  time. 
It  should  be  noted  that  with  the  application  of  the  A Tc  SH  correction,  the  variable  Tc  is  no  longer 
a  constant  over  the  globe  and  hence  can  no  longer  be  defined  as  the  nighttime  minimum  value 
of  Too.  The  choice  of  correcting  Tc,  as  opposed  to  T^  which  is  already  a  global  field,  is  the  con¬ 
vention  chosen  by  AFSPC,  and  one  to  which  we  adhere  in  order  to  directly  compare  our  model 
with  the  existing  model.  To  facilitate  comparison  of  the  individual  functions  of  an  SH  expansion 
with  those  of  other  functional  expansions  whose  indicies  increase  linearly,  we  define  the  ‘mode 
number’  to  be  used  throughout  this  chapter.  In  the  case  of  the  SH  functions,  each  mode  corre¬ 
sponds  to  a  single  term  within  Equation  fflin  the  order  of  the  summations,  with  the  cosine  terms 
preceeding  the  sine  terms  for  each  degree  and  order.  For  instance,  the  expansion  coefficients 
of  modes  1-4  of  the  SH  expansion  functions  would  be  C$,  C?,  C\ ,  and  S\,  respectively. 

Due  to  the  sparse  data  sets  available  for  near  real-time  estimation  of  the  SH  expansion  co¬ 
efficients,  it  was  necessary  to  truncate  the  SH  functions  in  Equation  [4]  for  degree  and  orders 
higher  than  [2x2]  (mode  9)  for  A Tc  and  [1  x  1  ]  (mode  4)  for  A Tx.  This  extension  has  come  to  be 
known  as  the  High  Accuracy  Satellite  Drag  Model  (HASDM).  The  goal  of  our  study  is  to  replace 
the  truncated  set  of  SH  functions  used  to  specify  A Tc  and  A Tx  with  a  global  basis  set  that  more 
accurately  represents  the  true  variability  of  the  thermosphere. 

2.2.2  The  TIE-GCM  First-Principles  Model 

The  Thermosphere-lonosphere-Electrodynamics  General  Circulation  Model  (TIE-GCM)  [59]. 
developed  by  the  National  Center  for  Atmospheric  Research  (NCAR),  is  a  comprehensive,  first- 
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principles,  three-dimensional,  non-linear  representation  of  the  coupled  thermosphere  and  iono¬ 
sphere  system  that  includes  a  self-consistent  solution  of  the  middle  and  low-latitude  dynamo 
field.  The  model  solves  the  three-dimensional  momentum,  energy  and  continuity  equations  for 
neutral  and  ion  species  at  each  time  step,  using  a  semi-implicit,  fourth-order,  centered  finite  dif¬ 
ference  scheme  on  each  pressure  surface  in  a  staggered  vertical  grid.  It  can  be  run  in  either 
serial  or  parallel  mode  on  a  variety  of  platforms,  including  Linux  workstations  and  supercomput¬ 
ers  [http://www.hao.ucar.edu/modeling/tgcm/tie.php,  accessed  24  Aug  2015]. 

The  default  model  time  step  of  120  seconds  is  used  when  solving  the  time-dependent  equa¬ 
tions;  however,  in  order  to  maintain  a  manageable  amount  of  data  the  output  is  only  recorded 
every  hour  for  the  purposes  of  this  study.  The  horizontal  grid  spacing  is  5°  in  latitude  and  lon¬ 
gitude,  with  a  vertical  spacing  of  a  half  scale  height.  In  this  configuration,  the  model  has  29 
pressure  surfaces  covering  the  altitude  range  from  97  km  to  500-700  km  depending  on  solar 
activity.  The  external  forcing  of  the  TIE-GCM  are  solar  irradiance  [166],  auroral  electron  precipita¬ 
tion  [61  ],  ionospheric  convection  driven  by  the  magnetosphere-ionosphere  current  system  [76], 
the  amplitudes  and  phases  of  tides  from  the  lower  atmosphere  pZI,  and  the  lower  boundary 
parameterization  of  eddy  diffusivity  [50]. 

2.2.3  Construction  of  Basis  Functions 

The  end  product  of  this  study  is  a  subroutine  that  operates  in  much  the  same  way  as  the 
J70  Dynamic  Calibration  Atmosphere  (J70DCA)  model  [[69].  However,  we  establish  an  improved 
set  of  basis  functions  to  more  efficiently  and  accurately  correct  the  J70  model  in  the  presence  of 
recent  satellite  data.  The  desired  basis  functions  should  have  the  ability  to  represent  the  behavior 
of  the  thermosphere  over  a  wide  range  of  geophysical  conditions.  To  that  end,  we  have  chosen 
to  simulate  such  conditions  using  the  TIE-GCM,  then  to  fit  the  output  to  the  low-dimensional  J70 
model,  and  finally  to  analyze  the  principal  components  of  the  fit  parameters.  The  remainder  of 
this  section  gives  an  in-depth  explanation  of  this  technique  and  a  description  of  the  resulting 
basis  functions. 

We  initially  ran  the  TIE-GCM  during  the  entire  year  of  2004.  However,  because  we  desire  a 
more  general  set  of  basis  functions,  and  not  one  capable  only  of  replicating  the  exact  thermo¬ 
spheric  conditions  present  in  a  free  run  of  the  TIE-GCM  during  a  single  year,  we  ran  an  ensemble 
of  the  TIE-GCM.  For  each  of  the  30  ensemble  members,  a  distinct  multiplicative  factor  is  applied 
to  the  F io  7  solar  flux,  the  Joule  heating,  and  the  eddy  diffusivity.  These  parameters  are  chosen 
to  account  for  any  natural  variability  of  the  input  forcing  and  lower  boundary  conditions,  as  well 
as  to  mitigate  any  internal  errors  within  the  model.  Within  each  ensemble  member,  a  random 
multiplicative  constant  is  applied  over  all  of  2004  for  each  of  the  three  adjustable  parameters. 
Because  the  multiplicative  factors  are  static  in  nature,  they  are  not  able  to  simulate,  for  instance, 
the  transient  response  of  the  thermosphere  during  geomagnetically  disturbed  times.  Instead, 
this  behavior  is  reproduced  by  the  input  timeseries  of  Interplanetary  Magnetic  Field  (IMF)  obser¬ 
vations  coupled  with  the  high-latitude  ionospheric  convection  model  of  Weimer  [I76[|.  The  role  of 
the  adjustable  Joule  heating  parameter,  however,  is  to  enable  the  ensemble  to  simulate  storms 
both  smaller  and  larger  than  those  that  actually  occurred  during  2004. 

While  the  nature  of  solar  flux  variability  can  be  quantified  given  the  relatively  long  record  of 
historical  data,  variability  of  the  other  two  chosen  parameters  is  somewhat  unobserved.  Because 
of  this,  the  choice  of  a  probability  density  function  (PDF)  from  which  to  draw  these  random  en- 
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Figure  1:  Distribution  of  ensemble  parameters.  Ensemble  input  parameters  F10.7  (top),  Joule 
heating  (middle),  and  eddy  diffusivity  (bottom)  after  application  of  randomly  generated  multi¬ 
plicative  factors.  The  histogram  (top  only)  shows  historical  data  from  2004-like  years  within  the 
last  5  solar  cycles,  black  lines  show  the  corresponding  PDFs,  and  color-coded  dots  indicate  the 
random  values  generated  for  each  ensemble  member.  The  histogram  and  all  PDFs  have  been 
normalized  such  that  the  integrated  area  under  the  curve  is  equal  to  1 . 


semble  parameters  is  somewhat  arbitrary.  Nevertheless,  certain  properties  were  desired  when 
choosing  the  PDF  of  the  multiplicative  factors.  Among  these  were:  (1 )  coverage  over  the  interval 
(0,oo),  (2)  the  ability  to  represent  skewness  in  data,  and  (3)  probability  Pr[X>a]  =  Pr[X<l/a] 
for  all  a>l.  The  second  property  gives  the  ability  to  more  closely  represent  the  positive  skew¬ 
ness  seen  in  historical  solar  flux  data  during  a  given  phase  of  the  solar  cycle.  The  third  property 
ensures,  for  example,  that  the  probability  of  generating  a  multiplicative  value  of  1/2  or  smaller 
is  the  same  as  that  of  generating  a  value  of  2  or  larger.  Given  these  desirable  properties,  a 
natural  choice  for  the  PDF  is  a  zero-mean  log-normal  distribution:  lnA/”(/i=0,  a2).  In  the  case  of 
solar  flux,  the  log-normal  parameter  a  can  be  fit  to  the  timeseries  of  F10.7  from  the  past  5  solar 
cycles  during  years  which  correspond  to  the  same  phase  of  the  solar  cycle  as  2004.  Using  the 
maximum  likelihood  technique,  a  standard  deviation  of  cr=0.20  was  estimated.  The  top  panel  of 
Figure [Bshows  the  histogram  of  historical  data  for  2004-like  years  along  with  the  resulting  PDF 
of  simulated  F10.7  calculated  for  a  day  in  which  the  measured  solar  flux  is  equal  to  the  median 
of  the  selected  historical  data  (i.e.  F10  7=1 18.8  sfu).  The  log-normal  PDF  fits  the  historical  data 
reasonably  well,  accounting  for  the  extended  upper-tail  as  well  as  the  sudden  drop-off  near  80 
sfu.  The  color-coded  dots  show  the  values  of  F107  that  would  be  generated  for  the  ensemble 
during  this  day. 
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Similar  arguments  were  used  in  the  construction  of  the  remaining  two  PDFs.  However,  with¬ 
out  the  aid  of  data,  there  was  little  hope  of  establishing  an  accurate  PDF  for  the  Joule  heating 
or  eddy  diffusivity  muliplicative  factors.  Therefore,  the  same  a  parameter  was  used  for  all  three 
PDFs.  In  theory,  the  a  parameters  for  the  latter  two  PDFs  could  be  adjusted  to  yield  the  best 
agreement  with  neutral  density  or  satellite  drag  data;  however,  this  process  was  deemed  com¬ 
putationally  impractical.  Moveover,  the  exact  shape  of  the  PDF  should  not  have  a  critical  effect 
on  the  overall  structure  of  the  resulting  basis  functions.  The  middle  panel  of  Figure  3]shows  the 
PDF  and  randomly  generated  ensemble  values  for  the  Joule  heating  multiplicative  factor,  taking 
into  account  the  TIE-GCM’s  default  value  of  1.5.  Finally,  the  lower  panel  shows  the  same  data 
for  eddy  diffusivity  during  a  day  in  which  the  default  value  equals  the  baseline  of  the  model  of 
Qian  et  al.  [S|  (i.e.  KE=A. 06xl0_6s_1). 

After  simulating  a  wide  range  of  conditions,  we  then  fit  the  TIE-GCM  output  neutral  mass 
density  to  the  J70  model  making  use  of  the  adjustable  parameters  A Tc  and  A Tx.  To  accomplish 
this,  J70  is  used  in  a  similar  fashion  to  the  method  of  Storz  et  al.  [69]:  however,  instead  of 
estimating  SH  expansion  coefficients,  A Tc  and  A Tx  are  estimated  directly  at  each  latitude/local 
time  grid  point  of  the  TIE-GCM  using  the  vertical  profile  of  mass  density  from  200  km  to  the 
upper  boundary  of  the  TIE-GCM.  Here  the  native  longitude  grid  of  the  TIE-GCM  is  transformed 
into  a  local  time  grid  using  the  simple  equation  6  =  UT  +  lon/15  (mod  24),  ignoring  the  +1 6/-1 4 
minute  offset  of  local  noon  over  the  course  of  the  year  [65].  (Note:  this  omission  can  easily 
be  remedied  in  future  model  versions  by  performing  this  fitting  procedure  in  the  appropriate 
functional  domain  as  opposed  to  the  gridded  domain).  Figure  [2]  shows  this  basic  process  for 
one  time  step  within  one  ensemble  member  and  at  one  location  in  latitude  and  local  time.  An 
iterative  nonlinear  least-squares  fit  is  employed  to  bring  the  original  J70  log-density  profile  into 
agreement  with  the  TIE-GCM  log-density  vertical  profile.  This  procedure  requires  knowledge  of 
the  partial  derivatives  dp/dTc  and  d p/dTx  with  respect  to  the  J70  model  at  each  location  and 
height,  which  are  calculated  by  the  J70DCA  algorithm  [69].  For  each  TIE-GCM  model  time  step 
and  ensemble  member,  this  produces  a  map  of  A Tc  and  A Tx  on  the  latitude/local  time  grid  of 
the  TIE-GCM  (such  as  those  shown  in  the  left  panels  of  figures  [6] and  [7)- 

It  should  be  noted  that  the  fits  of  the  TIE-GCM  vertical  profile  to  J70  are  not  perfect.  In  par¬ 
ticular,  there  is  no  guarantee  that  improving  the  agreement  of  the  vertical  mass  density  profiles 
between  J70  and  TIE-GCM  will  also  improve  the  agreement  of  the  vertical  temperature  pro¬ 
files.  This  behavior  is  caused  by  several  of  the  assumptions  of  J70  that  are  inconsistent  with 
the  TIE-GCM.  Among  these  is  J70’s  governing  assumption  of  diffusive  equilibrium  above  105 
km,  which  is  not  strictly  obeyed  by  the  TIE-GCM.  Photochemical  production/loss,  turbulent  diffu¬ 
sion  in  the  lower  thermosphere,  upwelling  [60]  and  horizontal  transport  [B  0]  all  commonly  cause 
departures  from  diffusive  equilibrium  within  the  TIE-GCM.  Additional  J70  assumptions  of  conse¬ 
quence  include  a  constant  lower  boundary  condition  for  mass  density,  temperature  and  mixing 
ratios  at  90  km,  an  imposed  static  model  describing  the  dissociation  of  molecular  oxygen  below 
105  km,  an  instantaneous  transition  from  turbulent  mixing  to  diffusive  equilibrium  at  105  km,  and 
a  temperature  profile  that  obeys  Equation  [3]  Therefore,  in  the  context  of  this  study,  the  adjusted 
temperature  profile  of  J70  should  be  regarded  as  a  mathematical  construct  whose  physical  inter¬ 
pretation  is  limited  by  the  degree  to  which  the  TIE-GCM  mimics  the  aforementioned  assumptions 
at  any  given  time. 

Through  analysis  of  the  A Tc  and  A Tx  corrections,  we  can  capture  the  most  significant  modes 
of  variability  using  the  principal  component  analysis  technique  [see  49].  Principal  component 
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Figure  2:  Temperature/density  vertical  profil  .  Temperature  (left)  and  density  (right)  vertical 
profiles  for  the  TIE-GCM  (black  squares,  density  only),  J70-unadjusted  (blue),  and  J70-adjusted 
(red)  at  one  latitude/local  time  grid  point.  Circles  indicate  unadjusted  and  adjusted  Tx  at  125  km 
and  ^  near  the  upper  boundary. 


analysis  is  a  linear  technique  that  provides  the  optimal  set  of  basis  functions  with  which  to  rep¬ 
resent  the  covariance  matrix  of  a  data  set  or  model  output.  Rearranging  the  A Tc  and  A Tx  fields 
into  an  m-by-n  matrix  called  F,  m  being  the  number  of  time  steps  (m=366  daysx24  hoursx30 
ensembles=263,520)  and  n  being  the  number  of  unique  grid  point  locations  (n=36  Iatitudesx72 
local  times=2,592),  we  compute  the  n-by-n  quasi-covariance  matrix  R  =  FTF.  At  this  point, 
we  solve  the  matrix  eigenvector/eigenvalue  problem:  Rfy  =  ^A,  where  ^  is  a  matrix  containing 
the  column  eigenvectors,  of  R,  and  A  is  a  diagonal  matrix  formed  from  the  eigenvalues,  A*, 
of  R.  The  eigenvectors,  are  referred  to  as  Principal  Component  (PC)  expansion  functions, 
while  the  eigenvalues  relate  to  the  amount  of  variance  of  the  original  fields  that  is  captured  by 
each  corresponding  eigenvector.  This  amounts  to  finding  the  set  of  orthonormal  basis  functions 
that  maximizes  the  projection  of  the  row-vectors  of  F  onto  each  basis  function.  In  addition  to 
providing  an  optimal  represention  of  the  covariance  matrix,  R,  a  coefficient  expansion  with  all  n 
eigenvectors  is  also  capable  of  representing  any  data  on  the  latitude/local  time  grid,  including  all 
data  contained  in  the  F  matrix.  This  is  simply  a  direct  consequence  of  the  orthonormality  of  the 
eigenvectors  and  could  be  said  of  any  basis  set  defined  on  the  same  grid.  However,  the  useful¬ 
ness  of  the  eigenvectors  comes  from  the  fact  that  most  of  the  variability  in  F  can  be  recovered 
using  only  a  handful  of  the  eigenvectors,  which  cannot  be  said  of  many  other  basis  sets. 
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2.3  Results  and  Discussion 


The  PC  procedure  is  carried  out  separately  for  the  A Tc  and  A Tx  fields.  Figures  [3]  and  8J 
show  the  first  9  PC  expansion  functions  for  A Tc  and  the  first  4  PC  expansion  functions  for  ATX, 
respectively,  after  rearranging  into  the  latitude/local  time  coordinate  frame.  The  different  modes 
are  a  combination  of  the  true  thermospheric  variability  with  the  error  in  both  the  TIE-GCM  and 
J70.  Much  of  the  variability  caused  by  physical  processes  that  cannot  be  captured  by  the  original 
J70  model  -  even  when  corrected  by  a  truncated  SH  expansion  -  can  be  represented  by  the  PC 
expansion  functions.  More  importantly,  the  cumulative  fraction  of  the  variance,  R,  captured  by 
each  mode  approaches  1  much  more  quickly  for  the  PC  expansion  than  it  does  for  the  SHs,  as 
shown  by  figure  [5}  In  fact,  if  we  were  to  attempt  to  capture  the  same  amount  of  variance  as 
the  first  PC  expansion  function  of  A Tc  and  A Tx  using  a  truncated  SH  expansion,  we  would  re¬ 
quire  the  expansion  functions  through  degree  and  order  [6x0]  (37  modes)  and  [2x0]  (5  modes), 
respectively. 

The  first  mode  of  PC  expansion  functions  accounts  for  98.1  and  87.8%  of  the  total  variance  in 
the  TIE-GCM  A Tc  and  A Tx  maps,  respectively,  as  shown  in  Figure  [51  Both  expansion  functions 
have  a  strong  mean  bias,  indicative  of  the  ability  to  correct  the  J70  model  bias.  In  general,  the 
PC  expansion  functions  are  a  convolution  of  model  error  -  introduced  by  the  TIE-GCM  as  well  as 
by  J70  -  with  the  true  thermospheric  variability,  and  are  not  necessarily  physically  significant  in 
spite  of  their  mathematical  significance.  However,  some  information  on  thermospheric  variability 
can  still  be  gleaned  from  several  of  the  lowest  modes.  For  instance,  modes  1,  3  and  4  of  A Tc 
(Figure [3)  have  a  strong  diurnal  component  at  low  to  mid-latitudes.  Mode  2  has  a  strong  semid¬ 
iurnal  component  evident  at  low  to  middle  latitudes,  giving  way  to  a  latitudinally  antisymmetric 
seasonal  component  near  the  poles.  Modes  such  as  6  and  7  give  the  appearance  of  having 
a  strong  diurnal  variation  that  maximizes  near  60°  latitude  in  the  southern  hemisphere.  In  this 
region,  the  phase  offset  between  the  two  modes  is  approximately  6  hours  in  local  time,  which 
allows  for  any  phase  and  amplitude  to  be  represented  by  the  appropriate  linear  combination. 
Given  their  non-stationary  nature  with  respect  to  local  time,  these  functions  are  more  likely  to 
be  associated  with  the  longitude-UT  effect  [seeSl|  related  to  the  offset  between  the  geographic 
and  geomagnetic  poles,  an  effect  that  isn’t  well  represented  by  the  longitudinally-independent 
J70  model  formulation  or  by  HASDM.  For  A Tx  (FigureS),  mode  1  has  a  strong  diurnal  compo¬ 
nent,  however  the  latitudinal  structure  is  unique  from  modes  1,  3  or  4  of  the  A Tc  PC  expansion 
functions.  Mode  2  for  A Tx  is  strikingly  similar  to  mode  2  for  A Tc,  but  with  an  increased  latitudinal 
range  of  the  semidiurnal  feature  and  an  antisymmetric  feature  that  is  subdued  in  the  southern 
hemisphere. 

For  an  initial  comparison  of  the  performance  between  the  new  and  currently-used  basis  func¬ 
tions,  we  reconstruct  the  A Tc  and  A Tx  fields  from  one  of  the  TIE-GCM  ensemble  members  using 
both  methods.  In  order  to  show  improvement  in  neutral  mass  density  specification,  we  also  re¬ 
construct  the  density  at  an  altitude  of  400  km  using  both  methods.  The  [2x2]  SH  representation 
of  A Tc  uses  9  orthogonal  functions  while  the  [1  x  1  ]  SH  representation  of  A Tx  uses  4  orthogonal 
functions.  To  stay  consistent  with  the  approach  used  by  HASDM,  we  also  limit  the  PC  expansion 
to  9  functions  representing  A Tc  and  4  functions  representing  A Tx.  At  each  time  step  of  the  TIE- 
GCM  ensemble  member,  the  truncated  set  of  expansion  functions  can  be  fit  to  the  A Tc  and  A Tx 
fields  using  a  least  squares  method  to  estimate  the  expansion  coefficients,  Cu  in  the  following 
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Figure  3:  ATC  expansion  functions.  First  9  modes  of  the  PC  expansion  functions  of  A Tc 
estimated  for  the  TIE-GCM  during  2004.  Each  PC  expansion  function  is  normalized  such  that 

VfVi  =  I- 


9 

Approved  for  public  release;  distribution  is  unlimited. 


ATX  PC  1  ATX  PC  2 


0  4  8  12  16  20  24  0  4  8  12  16  20  24 

Local  Time  (hr)  Local  Time  (hr) 


Figure  4:  ATX  expansion  functions.  First  4  modes  of  the  PC  expansion  functions  of  A Tx 
estimated  for  the  TIE-GCM  durinq  2004.  Each  PC  expansion  function  is  normalized  such  that 

=  1- 


formula: 

©M)  =  Ec<e<(^0)  (5) 

i 

where  0  is  the  reconstructed  field  (representing  either  A Tc  or  A Tx),  the  0,’s  are  the  expansion 
functions  (either  PCs  or  SHs),  ip  is  latitude,  6  is  local  time,  and  the  subscript  i  refers  to  the 
mode  of  the  expansion  function.  Density  at  400  km  is  then  calculated  from  J70  modified  by  the 
reconstructed  A Tc  and  A Tx  fields.  Figure  [6]  shows  the  original  A Tc  field  estimated  using  the 
TIE-GCM  neutral  density  as  input  data,  as  well  as  the  fits  using  modes  1-9  of  the  PC  and  SH 
expansions.  Figure [Tjshows  a  similar  fit  for  ATX  using  modes  1-4.  Figure!] shows  the  total  mass 
density,  p,  at  400  km  as  simulated  by  the  TIE-GCM  and  reconstructed  using  the  PC  and  SH 
expansion  functions.  The  epoch  for  figures and!] was  chosen  such  that  the  performance  of 
the  truncated  SH  expansion  is  typical  of  all  2004. 

To  compare  the  PC  reconstructions  to  the  baseline  SH  reconstructions  of  A Tc,  A Tx  and  p 
at  400  km,  we  employ  a  relative  Root  Mean  Squared  error  (rRMSe)  metric,  expressed  as  a 
percentage  increase  or  decrease  in  the  RMS  error: 


rRMSe 


-100  x 


Yip  Ye  (©pc(vj,  9)  —  ©TIE-GCM (<A  6))2 

Yv  Ye  (@sh (<£,  0)  -  0tie-gcm(^,  9))2 


(6) 


where  0PC  and  0Sh  represent  the  reconstructed  fields  (either  A Tc,  A Tx  or  p  at  400  km)  using  the 
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Figure  5:  Cumulative  fractional  variance.  Cumulative  fractional  variance  of  the  A Tc  (solid) 
and  A Tx  (dashed)  TIE-GCM  fields  that  is  recovered  by  the  PC  (blue)  and  SH  (green)  expansion 
functions  by  mode  number.  Note  that  the  first  25  SH  modes  are  equivalent  to  an  expansion 
through  degree  and  order  [4x4], 
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Figure  6:  ATC  reconstruction.  A Tc  estimated  for  the  TIE-GCM  for  a  typical  time  epoch  in  2004 
(left),  the  reconstruction  of  A Tc  using  the  9  lowest  order  PC  expansion  functions  (center),  and 
the  reconstruction  of  A Tc  using  a  [2x2]  (9  modes)  SH  expansion  functions  (right). 
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Figure  7:  ATX  reconstruction.  A Tx  estimated  for  the  TIE-GCM  for  a  typical  time  epoch  in 
2004  (left),  the  reconstruction  of  A Tx  using  the  4  lowest  order  PC  expansion  functions  (center), 
and  the  reconstruction  of  A Tx  using  a  [1  xl]  (4  modes)  SH  expansion  functions  (right). 


TIEGCM  p@  400km  (kg/m3):  Day  81, 23:00UT,  2004  x10-12  PC  p  @  400km  (kg/m3):  Day  81, 23:00UT,  2004  X10-12  SH  p  @  400km  (kg/m3):  Day  81, 23:00UT,  2004  x  iq'12 
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Figure  8:  Neutral  Density  reconstruction.  Neutral  density,  p,  at  400  km  estimated  for  the 
TIE-GCM  for  a  typical  time  epoch  in  2004  (left),  and  the  reconstructions  of  p  at  400  km  using  PC 
expansion  functions  (center)  and  SH  expansion  functions  (right). 


subscripted  method,  and  ©tie-gcm  represents  the  corresponding  field  as  calculated  directly  from 
the  TIE-GCM.  The  negative  sign  is  used  so  that  positive  rRMSe  values  indicate  an  improvement 
over  the  performance  of  the  SH  expansions.  For  the  reconstructions  of  A Tc  in  figure  [61  the  RMS 
error  improved  by  52.5%  when  using  the  PC  approach  instead  of  the  SH  approach.  For  the 
reconstructions  of  A Tx  in  figure  0  the  RMS  error  improved  by  38.7%.  As  for  the  reconstructions 
of  p  at  400  km  altitude  shown  in  figure  ®  the  RMS  error  improved  by  49.0%  during  this  time 
period. 

Figure  M  shows  the  27-day  smoothed  rRMSe  throughout  2004  for  the  A Tc  and  A Tx  fields. 
On  average,  the  PC  reconstructions  of  A Tc  and  A Tx  reduce  the  RMS  by  44.8%  and  45.1%, 
respectively,  over  the  SH  reconstructions.  Also  shown  in  Figure [9]is  the  27-day  smoothed  rRMSe 
for  p  at  400  km.  The  PC  reconstructions  of  p  reduce  the  RMS  by  32.9%  on  average  for  all  of 
2004. 
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PC  RMS  Error  Improvement  Relative  to  SH  expansions 
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Figure  9:  Relative  RMS  error.  27-day  smoothed  RMS  error  of  the  PC  reconstructions  of  p  at 
400  km  (black  line),  A Tc  (grey  line),  and  A Tx  (grey  dashed  line)  relative  to  RMS  error  of  the  SH 
reconstructions.  Positive  values  indicate  that  the  PC  reconstructions  are  outperforming  the  SH 
reconstructions. 


2.4  Conclusions 

We  have  presented  a  new  set  of  basis  functions  capable  of  representing  the  most  important 
modes  of  thermospheric  variability  within  the  framework  of  an  adjusted  static  diffusion  model. 
Our  approach  demonstrates  the  potential  for  significant  upgrades  to  current  operational  satellite 
drag  modeling  capabilities,  while  circumventing  many  of  the  technical  and  practical  issues  cur¬ 
rently  limiting  data  assimilation  with  physics-based  models.  Thermospheric  variability  is  specified 
by  the  TIE-GCM,  and  thus,  several  limiting  assumptions  should  be  pointed  out.  First,  the  simpli¬ 
fying  assumptions  of  hydrostatic  equilibrium  and  constant  gravity  are  imposed.  Secondly,  several 
physical  processes  are  not  fully  accounted  for,  such  as  the  influences  of  the  lighter  neutral  and 
ion  species,  [H]  and  [He].  In  addition,  eddy  diffusivity  at  the  lower  boundary  is  specified  by  an 
empirical  parameterization.  In  spite  of  these  and  other  simplifying  assumptions  along  with  the 
subsequent  error  that  they  impose  on  the  TIE-GCM,  the  modes  of  variability  of  the  TIE-GCM 
are  more  realistic  than  any  existing  empirical  model.  The  purpose  of  this  study  is  to  extract  the 
most  important  of  these  modes  and  use  them  to  efficiently  calibrate  empirical  models,  without 
the  increased  data  coverage  and  overhead  that  would  be  required  to  calibrate  a  high-dimension 
general  circulation  model  such  as  the  TIE-GCM. 

The  validation  presented  in  this  chapter  is  only  a  first  step.  A  comparison  of  the  new  and 
currently-used  techniques  using  actual  satellite  tracking  data  will  be  required  before  this  basis 
set  can  be  considered  validated.  As  mentioned  in  Section  3,  several  of  the  PC  modes  exhibit 
traits  which  are  not  included  in  the  J70  model,  e.g.  the  longitude-UT  effect.  Thus  we  are  able 
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to  include  and  correct  for  modes  of  thermospheric  variability  not  captured  the  current  implemen¬ 
tation  of  HASDM.  However,  some  of  these  traits  may  not  be  observable  from  the  ground-based 
satellite  tracking  data  set  currently  used  to  drive  HASDM.  Validation  studies  with  these  ground- 
based  satellite  tracking  data  is  needed  to  resolve  such  issues,  as  well  as  guide  any  necessary 
modifications  to  the  PCs  to  provide  the  needed  upgrade  to  operational  satellite  drag  modeling 
capabilities. 
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3  First-Principles  Helium  Model 


3.1  Background 

The  presence  of  helium  as  a  major  component  in  the  Earth’s  upper  thermosphere  and  lower 
exosphere  was  first  inferred  from  measurements  of  satellite  drag.  By  analyzing  orbital  variations 
of  the  Echo  1  satellite  orbiting  above  1000  km,  Nicolet  [[47[]  reasoned  that  atomic  oxygen  was 
incapable  of  producing  the  observed  satellite  deceleration  given  reasonable  values  of  exospheric 
temperature.  Likewise,  atomic  hydrogen  concentrations  were  thought  to  be  much  too  low  to 
create  such  a  deceleration. 

Increasingly  direct  evidence  of  helium’s  presence  soon  emerged  from  in  situ  mass  spectrom¬ 
eter  measurements  taken  onboard  Explorer  17  [54].  Concomitant  with  this  confirmation  was  the 
hint  of  a  significant  seasonal-latitudinal  variation  in  the  helium  distribution,  relative  to  the  other 
measured  constituents  (i.e.  molecular  nitrogen  and  atomic  oxygen).  Soon  thereafter,  strong 
semi-annual  variations  inferred  from  the  satellite  drag  acting  on  Echo  2  [6]  around  1 1 00  km  were 
linked  to  seasonal  variations  of  helium  concentration.  Keating  and  Prior  [26]  confirmed  this  result 
with  satellite  drag  data  from  the  Explorer  9,  19,  and  24  satellites.  They  also  noted  an  apparent 
enhancement  near  the  winter  pole,  which  they  termed  the  “winter  helium  bulge,”  with  an  approx¬ 
imate  winter-to-summer  ratio  of  2.5.  Subsequent  drag-inferred  calculations  by  Keating  etai.  [27] 
yielded  ratios  in  excess  of  3  at  an  altitude  of  850  km. 

Reber  et  al.  [[55|,  using  mass  spectrometer  measurements  from  the  Ogo  6  satellite,  showed 
an  order-of-magnitude  difference  between  the  helium  content  in  winter  and  summer  hemispheres 
near  400-600  km  altitude.  This  disagreement  with  previous  results  highlighted  the  limitations  of 
the  drag-inferred  technique,  specifically,  reliance  on  the  assumption  of  diffusive  equilibrium  to 
separate  composition-induced  mass  density  variations  from  those  caused  by  temperature.  In 
response,  Keating  et  al.  [28]  augmented  their  drag-inferred  technique  to  include  a  description  of 
the  background  composition  that  was  consistent  with  the  available  mass  spectrometer  data.  New 
ratios  in  excess  of  an  order  of  magnitude  could  then  be  obtained  through  this  method  as  well. 
In  addition  to  establishing  a  larger  bulge  ratio,  Reber  et  al.  [t55^  noted  a  strong  correlation  of  the 
maximum  helium  density  with  the  location  of  the  winter  geomagnetic  pole.  This  was  interpreted 
as  a  sensitivity  of  the  helium  distribution  to  the  thermospheric  wind  system. 

In  addition  to  high-latitude  variations  near  the  solstices,  Newton  et  al.  [Si]  detected  a  strong 
local  time  preference  for  helium  concentration  as  measured  by  mass  spectrometers  on  the  low- 
inclination  San  Marco  3  satellite.  Reber  etai.  [156]  and  Mayretal.  |44]|  discussed  similar  variations 
manifest  within  the  Ogo  6  density  model  [16].  These  findings  showed  a  preference  of  the  diurnal 
maxima  toward  earlier  times  for  species  with  small  molecular  masses,  with  the  opposite  being 
true  for  species  of  large  mass.  The  San  Marco  3  observations,  taken  at  altitudes  near  225  km, 
showed  a  preference  toward  the  06-09  LT  sector  while  those  taken  by  Ogo  6,  near  450  km, 
showed  maxima  closer  to  10  LT. 

The  realization  of  these  phenomena  motivated  several  modeling  studies  to  uncover  the  mech¬ 
anism  responsible  for  the  counterintuitive  distribution  of  helium  in  the  thermosphere.  Noticing 
that  helium  vertical  profiles  measured  by  several  rocket-based  mass  spectrometers  departed 
quite  drastically  from  diffusive  equilibrium,  an  early  study  by  Kasprzak  [f25|  invoked  an  addi¬ 
tional  diffusive  flux  in  order  to  reconcile  the  observations  with  his  model.  This  treatment  required 
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vertical  fluxes  on  the  order  of  6  x  10s  and  2  x  1010  cm_2s_1  for  summer  and  winter  conditions, 
respectively,  over  an  altitude  range  of  120-200  km.  Kockarts  [31]  later  noted,  however,  that 
these  values  were  larger  than  the  maximum  flux  allowed  by  molecular  diffusion,  thus  requiring 
an  additional  mechanism  of  transport. 

Johnson  and  Gottlieb  [23]  used  basic  considerations  of  continuity  to  show  that  a  general 
summer-to-winter  flow  of  the  major  atmospheric  constituents  could  account  for  a  buildup  of  he¬ 
lium  in  the  winter  polar  regions.  Without  discounting  these  findings,  several  attempts  were  made 
to  ascertain  the  effect  of  atmospheric  fluctuations  on  helium  transport.  Hodges  [19]  modeled 
large-scale  fluctuations  as  monochromatic  plane  waves,  which  effected  a  downward  transport 
and  an  overall  decrease  to  the  scale  height  of  species  with  masses  smaller  than  the  mean 
mass.  Similarly,  Kockarts  [H]  derived  the  eddy  diffusivity  profile  necessary  to  reconstruct  the 
winter  helium  bulge  observations  of  Reber  et  al.  [tUf],  under  the  assumption  of  molecular  diffu¬ 
sion  in  the  absence  of  wind.  Results  from  these  studies  suggested  that  eddy  diffusion  could  in 
fact  control  the  global  helium  distribution.  However,  recreating  the  observed  winter  bulge  ratios 
required  more  than  an  order-of-magnitude  increase  in  eddy  diffusivity  from  winter  to  summer 
hemispheres.  These  results  were  qualitatively  consistent  with  each  other,  yet  they  implied  that 
similar  latitudinal  signatures  should  be  evident  in  other  minor  atmospheric  constituents,  a  feature 
that  was  inconsistent  with  previous  observations  of  atomic  oxygen  [31J. 

Reber  and  Hays  [53]  performed  a  more  rigorous  treatment  of  the  effects  of  circulation  on  the 
distribution  of  helium.  Included  in  their  model  were  the  effects  of  molecular  and  eddy  diffusion 
as  well  as  a  parameterized  circulation  pattern  of  the  background  gas  that  satisfied  continuity 
requirements  and  could  be  tuned  to  simulate  a  given  level  of  summer-to-winter  flow.  Combining 
the  equations  of  continuity  and  momentum  for  a  minor  species  led  to  an  accurate  representation 
of  previous  winter  helium  bulge  observations.  The  idea  that  the  winter  helium  bulge  could  be 
completely  explained  by  seasonal  circulation  patterns  led,  however,  to  an  apparent  paradox. 
At  times  of  high  solar  flux,  when  an  enhanced  summer-to-winter  flow  had  been  expected  to 
occur,  smaller  pole-to-pole  helium  ratios  had  been  observed.  Reber  and  Hays  [53]  explained  the 
discrepancy  by  invoking  the  mechanism  of  exospheric  flow,  whereby  during  times  of  high  solar 
flux,  increased  temperatures  in  the  upper  thermosphere  drive  a  larger  exospheric  flow  directed 
away  from  the  winter  bulge.  The  balance  between  the  circulation-induced  effects  and  exospheric 
transport  was  found  to  control  the  magnitude  of  the  latitudinal  gradient  in  helium  concentration 
that  could  be  supported  by  the  atmosphere. 

By  analyzing  the  combined  equations  of  continuity  and  momentum  for  a  minor  species,  Reber 
and  Hays  [53]  and  Hays  et  al.  [12]  identified  the  vertical  advection  term  as  being  responsible 
for  establishing  the  seasonal  distribution  of  helium.  In  the  presence  of  diffusively  separated 
atmospheric  constituents,  this  term  leads  to  increased  helium  densities  in  regions  of  downwelling 
and  decreased  densities  in  regions  of  upwelling.  The  opposite  behavior  is  implied  for  species, 
such  as  argon,  that  are  heavier  than  the  local  mean  mass.  Reber  [fll]  further  explained  that 
in  order  to  perturb  composition  from  the  distribution  prescribed  under  conditions  of  diffusive 
equilibrium,  the  vertical  winds  must  be  significant  in  relation  to  a  characteristic  vertical  diffusive 
velocity,  vD  =  D/H,  where  D  is  the  mutual  diffusion  coefficient  and  H  is  the  atmospheric  scale 
height. 

Contemporaneous  works  by  Mayr  and  Volland  [SI,  43]  asserted  a  similar  yet  distinct  per¬ 
spective  on  the  matter.  Mayr  et  al.  [45]  summarized  these  findings  and  coined  the  phrase 
“wind-induced  diffusion,”  describing  horizontal  transport  in  the  presence  of  diffusively  sepa- 
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rated  constituents.  Both  groups  agreed  that  the  interaction  between  helium  and  the  back¬ 
ground  circulation — consisting  of  upwelling  in  the  summer  hemisphere,  summer-to-winter  flow, 
and  downwelling  in  the  winter  hemisphere — would  lead  to  a  winter  helium  bulge  consistent  with 
observations.  However,  Reber  and  Hays  [53]  suggested  that  the  transport  mechanism  was  re¬ 
lated  to  the  vertical  advective  motion  in  the  presence  of  diffusive  separation,  while  Mayr  et  at. 
[45]  believed  horizontal  bulk  motion  in  the  presence  of  diffusive  separation  to  be  responsible. 

As  the  basic  mechanism  causing  the  observed  helium  behavior — i.e.  circulation  within  a 
diffusively  separated  atmosphere — continued  to  mature,  several  successful  satellite  mass  spec¬ 
trometer  missions  served  to  refine  these  theories  and  document  the  phenomenological  impli¬ 
cations.  The  open  source  mass  spectrometers  on  Atmospheric  Explorer  satellites  (AE-C  -D, 
and  -E)  were  used  by  several  investigators  to  further  quantify  seasonal  variations  [0,  4l],  §]. 
Reber  et  at.  [f57l  also  analyzed  these  data  to  study  waves  in  composition,  showing  coherent 
phase  relationships  between  the  various  constituents.  Hedin  and  Carignan  pj||  used  data  from 
the  Dynamics  Explorer  2  (DE-2)  satellite  to  show  that  even  during  geomagnetically  quiet  times, 
signatures  of  helium  depletion  are  present  near  the  magnetic  poles.  These  data  sets  now  com¬ 
prise  the  majority  of  our  understanding  of  thermospheric  composition,  the  empirical  basis  of 
which  is  embodied  by  the  Mass  Spectrometer  and  Incoherent  Scatter  (MSIS^series  of  models 


[13,  Ll4l,  148],  successors  of  the  Ogo  6  model.  More  recently,  Thayer  et  at.  [74]  inferred  strong 


signatures  of  helium  from  differences  in  total  mass  densities  measured  at  two  different  altitudes 
by  high-precision  accelerometers  on  board  the  Challenging  Mini-Satellite  Payload  (CHAMP)  and 
Gravity  Recovery  and  Climate  Experiment  (GRACE)  satellites  (7Q|].  Liu  et  at.  Q  extended  this 
work,  showing  that  the  response  of  the  mass  density  vertical  profile  during  a  geomagnetic  dis¬ 
turbance  is  quite  sensitive  to  the  atomic  oxygen/helium  transition  height. 

The  remainder  of  this  chapter  is  organized  as  follows.  Section 3)2] introduces  a  self-consistent 
method  for  calculating  helium  abundances  and  transport  by  modifying  an  existing  general  circula¬ 
tion  model  of  the  thermosphere.  Unlike  previous  formulations,  we  do  not  impose  the  assumption 
that  helium  remains  a  minor  species  throughout  the  model  domain,  which  can  have  deleterious 
effects  at  high  altitudes.  Section  [3)3]  highlights  the  salient  features  of  the  new  model,  including 
helium’s  role  in  determining  mean  mass,  total  mass  density,  pressure  level  height  and  winds. 


3.2  Methods,  Assumptions,  and  Procedures 

3.2.1  The  TIE-GCM  Baseline  Model 


The  model  developments  described  in  this  chapter  have  been  applied  to  the  National  Cen¬ 
ter  for  Atmospheric  Research  Thermosphere-lonosphere-Electrodynamics  General  Circulation 
model  (NCAR/TIE-GCM)  v.1.95  [63,  [59|,  and  are  slated  for  inclusion  in  the  next  TIE-GCM  and 
TIME-GCM  [62[|  model  versions.  The  TIE-GCM  is  a  first-principles  upper  atmospheric  general 
circulation  model  that  solves  the  Eulerian  continuity,  momentum,  energy,  and  composition  equa¬ 
tions  for  the  coupled  thermosphere-ionosphere  system.  The  vertical  coordinate  is  specified  by 
log-pressure  levels  in  half-scale  height  increments,  providing  coverage  in  altitude  of  approxi¬ 
mately  97  km  to  600  km,  the  latter  being  dependent  on  solar  activity. 

Tidal  forcing  at  the  lower  boundary  is  specified  by  the  Global  Scale  Wave  Model  (TT||.  Annual 
and  semi-annual  variations  in  sub-grid  turbulent  fluctuations  are  taken  into  account  by  applying 
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seasonal  variation  of  the  eddy  diffusivity  coefficient  at  the  lower  boundary  [5Q,  51].  Based  on 
measurements  from  the  Mauna  Loa  Observatory  [[291],  the  mixing  ratio  of  C02  imposed  at  the 
lower  boundary  was  set  to  364  ppmv  for  1 996,  increasing  linearly  by  1 .5  ppmv  per  year  thereafter. 

In  the  simulations  presented  throughout  this  chapter,  solar  irradiance  is  specified  in  a  manner 
consistent  with  Solomon  et  at.  [67].  The  M10.7  index  is  used  in  place  of  the  F10.7  solar  proxy  in 
an  effort  to  better  capture  solar  UV  and  EUV  irradiance  during  the  deep  solar  minimum  of  2008. 
The  M10.7  index  derives  from  the  magnesium  core-to-wing  (Mgll  c/w)  of  Viereck  et  ai  [75]  via  a 
linear  fit  to  the  F10.7  proxy  calculated  during  1978-2007  [67].  With  this  normalization,  M10.7  can 
be  used  in  place  of  F10.7  to  drive  the  EUVAC  proxy  model  [see  58],  73  66]. 

Magnetospheric  inputs  to  the  polar  regions  are  specified  by  an  applied  electric  potential  pat¬ 
tern  and  an  auroral  precipitation  oval.  The  Heelis  et  ai  |[l8j]  empirical  specification  of  magneto- 
spheric  potential  in  the  ionosphere,  which  is  parameterized  by  the  3-hour  geomagnetic  KP  index, 
is  the  standard  TIE-GCM  input  and  is  employed  for  the  simulations  presented  throughout  this 
chapter.  Auroral  precipitation  is  applied  as  described  by  Roble  and  Ridley  [6l|  based  on  the 
estimated  hemispheric  power  of  precipitating  electrons.  The  empirical  estimate  of  this  power  as 
it  depends  on  KP  has  been  increased  from  its  original  formulation  by  a  factor  of  ~2,  based  on 
results  from  the  Global  Ultraviolet  Imager  (GUVI)  on  the  TIMED  satellite  [78]. 

The  TIE-GCM  uses  the  method  outlined  by  Richmond  et  at  [59]  to  calculate  the  low-latitude 
ionospheric  electrodynamo  driven  by  conductances  and  neutral  dynamics.  The  calculated  elec¬ 
tric  potential  is  merged  with  the  externally  imposed  potential  within  each  polar  cap,  using  cross¬ 
over  boundaries  that  vary  dynamically  with  the  size  of  the  magnetospheric  potential  pattern. 
See  Solomon  et  at  Q,  section  2.3,  for  further  detail  concerning  the  high-latitude  inputs,  and 
Solomon  et  ai  [l67],  section  4,  for  a  discussion  of  model  uncertainties. 


3.2.2  Implementing  Helium  as  a  Major  Species 

The  equations  describing  the  transport  and  concentration  of  the  various  components  in  the 
upper  atmosphere  are  coupled  to  one  another  through  diffusive  and  chemical  processes.  When 
solving  for  the  concentration  of  a  minor  species  [63],  several  terms  in  the  fully  coupled  compo¬ 
sition  equation  are  assumed  to  be  small.  With  the  neglect  of  these  terms,  the  solution  of  the 
major  species  composition  becomes  dynamically  decoupled  from  that  of  the  minor  species  com¬ 
position,  leading  to  a  more  efficient  segmented  numerical  solution.  The  main  terms  that  must 
be  neglected  are  those  in  the  diffusion  matrix  describing  the  acceleration  experienced  by  any 
major  species  caused  by  collisions  with  the  minor  species  as  well  as  those  that  account  for  the 
effect  that  the  minor  species  has  on  the  mean  mass  and  scale  height  of  the  atmosphere.  It  is 
straightforward  to  show  that  the  effect  of  these  terms  is  small  when  the  mass  mixing  ratio  of 
the  minor  species  in  question  is  also  small.  Helium  as  a  minor  species  in  the  TIE-GCM  was 
recently  implemented  by  Liu  et  ai  [35].  While  this  approach  demonstrated  the  model’s  ability  to 
accumulate  helium  in  the  winter  hemisphere,  it  required  the  ad  hoc  inclusion  of  helium  into  the 
scale  height  calculation  in  order  to  avoid  unrealistically  high  values  during  long  simulations. 

As  a  simple  test,  an  empirical  approach  can  be  used  to  ascertain  whether  or  not  helium 
satisfies  the  assumptions  necessary  to  treat  it  as  a  minor  constituent.  We  employ  the  MSIS 
model  [48].  which  represents  helium  abundance  in  an  averaged  sense  as  observed  by  mass 
spectrometer  observations  spanning  several  decades.  However,  care  must  be  taken  when  con¬ 
verting  between  the  vertical  coordinate  systems  of  MSIS  and  TIE-GCM.  The  TIE-GCM  uses 
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log-pressure,  z  =  ln(p0/p),  as  its  vertical  coordinate,  where  p0  is  a  reference  pressure  set  to 
5  x  10~4g/(cm  ■  s 2).  In  order  to  obtain  a  reasonable  estimate  of  the  amount  of  helium  that  should 
be  present  within  the  vertical  domain  of  the  TIE-GCM,  it  is  necessary  to  compute  MSIS  densities 
with  respect  to  the  TIE-GCM’s  log-pressure  scale.  Using  the  ideal  gas  law,  we  directly  calculate 
the  log-pressure  level  from  the  number  densities  and  temperatures  specified  by  MSIS. 

As  molecular  diffusion  becomes  dominant  with  increasing  height,  a  neutral  species  of  com¬ 
paratively  small  mass  such  as  helium  will  increase  in  relative  concentration.  Due  to  the  interac¬ 
tion  between  global  circulation  and  molecular  diffusive  flow,  the  largest  values  tend  to  occur  at 
high  latitudes  in  the  winter  hemisphere  [e.g.  53,  45].  Figure  [TO] shows  that  under  these  conditions 
and  near  the  top  level  of  the  TIE-GCM  (i.e.  roughly  500-700  km,  depending  on  solar  flux),  helium 
mass  mixing  ratios  exceed  0.8  during  solar  maximum  conditions  and  0.9  during  solar  minimum 
conditions.  Had  we  instead  queried  MSIS  using  the  geometric  heights  calculated  by  TIE-GCM 
as  our  vertical  coordinate,  values  just  below  0.5  would  have  been  obtained.  As  will  be  shown 
in  Section  [3T31  this  discrepancy  stems  from  an  underestimation  of  the  geometric  height  in  the 
upper  thermosphere  by  the  original  TIE-GCM  code  due  to  the  neglect  of  helium.  In  either  case, 
empirical  evidence  suggests  that  helium  becomes  a  major  neutral  component — and  perhaps  the 
dominant  component — under  certain  conditions  within  the  spatial  domain  of  the  TIE-GCM.  In 
light  of  these  findings,  the  remainder  of  this  section  covers  the  expansion  of  the  major  neutral 
species  composition  equation  and  other  modeled  processes  from  a  3-constituent  description  p] 
to  a  4-constituent  description  in  order  to  account  for  the  significant  effects  of  helium. 

The  evolution  of  the  major  neutral  species  composition  can  be  expressed  using  the  following 
vector  equation  (see  the  Appendix  for  derivation  and  a  complete  definition  of  variables): 
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dt 
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-i  a_ 
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rTbo\0-25 


m,No  V  T  / 
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|[e-^)(|  +  lf)  *] 
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The  meanings  of  several  variables  have  been  modified  from  those  originally  intended  by  Dickin¬ 
son  et  al.  [|9].  is  now  the  vector  of  mass  mixing  ratios  for  02,  O,  and  He,  while  the  mass  mixing 
ratio  of  the  remaining  major  constituent  N2  is  specified  by  -0tv2  =  1  -  -0o2  ~i>o  ~  V’/fe-  Molecular 
and  thermal  diffusion  are  accounted  for  by  the  first  term  on  the  right  side  of  Eq.  (7),  eddy  diffusion 
by  the  second,  horizontal  and  vertical  advection  by  the  third,  and  chemical  sources  and  sinks  by 
the  fourth. 

L  is  a  diagonal  matrix  operator  with  elements: 


L  =  —  -  (i  -  —  -  ——  -  in) 

dz  \  fh  fh  dz  T  dz ) 

which  have  been  expanded  to  describe  thermal  diffusion,  a  phenomenon  which  becomes  impor¬ 
tant  for  species  such  as  helium  whose  masses  are  quite  different  from  the  mean  mass.  We  use 
a  simplified  formulation  of  thermal  diffusion  that  is  analogous  to  its  appearance  in  the  binary  dif¬ 
fusion  equations,  after  Colegrove  et  al.  [5].  In  this  treatment,  a  constant  value  of  aHe  =  -0.38  is 
used.  While  this  value  is  characteristic  of  small  concentrations  of  helium  diffusing  through  molec¬ 
ular  nitrogen,  this  assumption  is  reasonably  accurate  at  altitudes  where  significant  temperature 
gradients  exist  (i.e.  below  ~200  km)  p]]. 
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Figure  10:  Mass  mixing  ratio  vertical  profil  .  MSIS  mass  mixing  ratios  for  02  (blue),  O 
(green),  N2  (red),  and  He  (cyan)  calculated  on  the  vertical  log-pressure  scale  in  the  vicinity  of 
the  winter  helium  bulge  for  solar  maximum  (solid  lines/black  altitude  labels,  21  Dec.,  2000)  and 
minimum  (dashed  lines/grey  altitude  labels,  21  Dec.,  2008)  conditions. 


The  normalized  molecular  diffusion  matrix,  a,  couples  the  major  components  to  one  another. 
As  can  be  seen  in  Eqs.  (34)  and  (39)  in  the  appendix,  the  strength  of  this  coupling  depends 
on  the  mutual  diffusion  coefficients.  Dickinson  et  al.  [9f|  assumed  these  coefficients  to  take  the 
form  D  =  L>o(T/Too)175(poo/p)  for  the  major  species,  after  Colegrove  et  al.  [[5j.  Accordingly, 
the  elements  of  a  have  been  normalized  by  this  functional  form.  Mutual  diffusion  coefficients 
between  helium  and  the  other  three  major  species  take  a  similar  form,  yet  with  exponents,  s,  that 
deviate  slightly  from  1 .75,  as  seen  in  Table  (3.2.2).  These  differences  have  been  accounted  for  by 
applying  correction  factors  of  the  form  (T/T00)L75”S  to  the  appropriate  terms  within  the  diffusion 
matrix  a.  In  the  absence  of  these  corrections,  the  coefficient  describing  the  mutual  diffusion 
between  helium  and  atomic  oxygen  would  remain  reasonably  accurate,  yet  those  describing  the 
interaction  of  helium  with  molecular  species  would  be  approximately  5%  low. 

The  chemical  source  and  sink  matrix,  s,  also  serves  to  couple  the  major  species  to  one 
another.  In  the  case  of  helium,  however,  all  chemical  and  photochemical  rates  have  been  set  to 
zero,  consistent  with  our  assumption  of  inertness.  Therefore,  our  current  model  implementation 
is  appropriate  for  the  study  of  the  dynamical  behavior  of  helium  as  an  ideal  inert  tracer. 

The  neutral  thermodynamic  properties  of  specific  heat,  cp,  molecular  viscosity,  km,  and  con¬ 
ductivity,  kt,  have  been  augmented  to  include  the  effects  of  helium.  The  following  equations  are 
now  used  [1J]: 
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Table  1:  Mutual  diffusion  coefficients  for  helium  with  major  species1. 


i-g 

a 

s 

He-02 

0.649 

1.710 

He-0 

0.866 

1.749 

He-N2 

0.622 

1.718 

lDy  —  a  (T/T00] 

f  (poo/p),  Too=273  K,  _p0o=106  g/(cm  ■ 

•  s2)  [cf.  1 ,  table  15.1]. 
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(56. 0(no2  +  n^)  +  75.9 no  +  299.0 nHe)  erg  ■  cm  1K  1 


n 


(9) 

(10) 

(11) 


where  R  is  the  universal  gas  constant,  T  is  the  neutral  temperature  in  units  of  Kelvin,  n*  refers  to 
the  number  density  of  the  subscripted  species,  and  n  is  the  total  number  density. 

Additionally,  in  the  description  of  ambipolar  diffusion,  the  collision  frequency,  uin,  has  been 
updated  to  account  for  nonresonant  collisions  between  0+  ions  and  neutral  helium  atoms.  The 
following  form  is  adopted  [64]: 


uin  =  1  x  10  10(6. 64 rto2  +  0.367 n0\/Tr{  1  —  0.064 log10 Tr )  +  6.82 +  1.32 nue)  (12) 

where  Tr  =  (Tt  +  T)/ 2  is  the  average  of  the  ion  and  neutral  temperatures.  Tr,  uin  and  n,  are  in 
units  of  Kelvin,  s-1  and  cm-3,  respectively. 


3.2.3  Boundary  Conditions 

At  the  lower  boundary  of  the  model,  atomic  and  molecular  oxygen  adhere  to  the  conditions 
specified  in  the  original  TIE-GCM  implementation,  namely,  that  the  peak  of  the  atomic  oxygen 
density  profile  lies  at  the  lower  boundary  and  the  total  amount  of  oxygen  atoms  remains  constant 
making  up  23.4%  of  the  total  mass.  In  addition,  we  specify  a  constant  lower  boundary  mass 
mixing  ratio  for  helium  of  1.154 x  10-6.  In  terms  of  mass  mixing  ratios,  these  considerations  take 
the  following  form:  (1)  di>0/dz  =  ip0,  and  (2)  ipo2  +^o  =  0.234,  and  (3)  tpHe  =  1.154 xlO-6. 

Near  the  upper  boundary  of  the  model,  either  atomic  oxygen  or  helium  typically  dominates 
the  composition,  depending  on  season,  solar  flux,  and  location.  The  original  upper  boundary 
of  the  TIE-GCM  is  specified  by  diffusive  equilibrium  for  neutral  species,  i.e.  LT  =  0.  However, 
the  large  thermal  velocity  of  helium  warrants  proper  consideration  of  helium  transport  processes 
occurring  above  the  upper  boundary  in  a  near-collisionless  environment.  While  the  classical 
thermal  escape  flux  of  helium  is  several  orders  of  magnitude  too  low  to  have  a  noticeable  effect 
on  the  global  helium  content,  the  lateral  transport  of  helium  atoms  with  ballistic  trajectories  is 
significant.  Hodges  and  Johnson  [2lJ  and  Hodges  d20t|  outline  a  method  for  approximating  this 
type  of  transport,  expressing  it  as  a  vertical  outward  particle  flux: 
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$  =  -V2  (nvH2P)  (13) 

where  V2  is  the  surface  Laplacian.  The  variables  $,  n,  v,  and  H  are  respectively  the  vertical 
particle  flux,  number  density,  mean  thermal  speed,  and  scale  height,  all  specific  to  helium.  P,  a 
dimensionless  factor  arising  from  integration  over  Maxwellian  distributions,  has  a  weak  depen¬ 
dence  on  neutral  temperature  that  can  be  adequately  approximated  by  [21]: 


+  (14) 

V  3300 )  v  ' 

for  neutral  temperature,  T,  in  units  of  Kelvin.  Inherent  in  these  equations  is  the  assumption  that 
collisions  do  not  occur  above  the  upper  boundary  of  the  TIE-GCM. 

In  practice,  this  vertical  flux  can  be  prescribed  at  the  upper  boundary  of  the  model  as  a 
diffusive  flow.  The  following  vector  equation  describing  molecular  diffusion  is  used: 


W  d  =  t 
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Too 
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0.25 


-^a-1  L* 


gmN2 


(15) 


where  wD  is  the  (3x1)  vector  of  vertical  diffusive  mass  flow  rates  for  02,  O,  and  He,  respectively. 

it  follows  that  the  diffusive  mass 


From  the  derivation  of  Eq.  05)  in  the  appendix  (see  Eq. 
fluxes  of  all  neutral  species  sum  to  zero.  Because  molecular  oxygen  and  nitrogen  densities  are 
small  near  the  upper  boundary,  we  enforce  this  constraint  by  assuming  that  any  outward  (inward) 
mass  flux  of  helium  is  balanced  by  an  inward  (outward)  flux  of  atomic  oxygen.  Any  error  that  this 
assumption  incurs  in  the  solution  of  atomic  oxygen  concentration  is  diminished  by  the  factor  of  4 
difference  between  the  mass  of  oxygen  and  helium  atoms. 

In  the  current  implementation  of  our  model,  the  argument  of  the  Laplacian  from  Eq.  03) 
is  transformed  into  a  non-aliasing  spherical  harmonic  expansion.  This  is  completed  using  the 
technique  of  Swarztrauber  [73],  modified  to  accommodate  the  TIE-GCM’s  horizontal  grid  which 
is  offset  from  the  pole  by  a  half-grid  increment.  The  flux,  <E>,  is  then  calculated  using  the  well- 
known  eigenfunction/eigenvalue  relation: 


V2ym  =  _n(n  +  l)c  (16) 

where  Y™  refers  to  the  spherical  harmonic  function  of  degree  n  and  order  m,  and  R  is  a  charac¬ 
teristic  radius  of  the  exobase.  In  the  current  implementation,  R  has  been  set  to  the  radius  of  the 
Earth  for  consistency  with  calculations  of  other  horizontal  derivatives  within  the  TIE-GCM.  The 
mass  flux  required  by  the  left-hand-side  of  Eq.  05)  can  then  be  obtained  by  transforming  back 
to  the  spherical  grid  and  multiplying  the  obtained  particle  flux  by  the  molecular  mass  of  helium. 
The  advantage  of  using  this  technique  in  place  of  finite  differences  for  calculating  the  Laplacian  is 
that  waves  are  resolved  uniformly  over  the  Earth.  Therefore,  the  growth  of  numerical  instabilities 
can  be  controlled  by  truncating  the  expansion  prior  to  transforming  back  to  the  spherical  grid. 
We  note  that  the  degree  of  truncation  required  is  sensitive  to  the  level  of  the  upper  boundary, 
the  grid-size,  and  the  time  step.  When  using  the  default  5°x  5°xH/2  spatial  grid  with  upper 
boundary  of  z  =  +7  and  a  120  second  time  step,  we  have  found  that  a  triangular  truncation  of 
degrees  higher  than  4  is  sufficient  to  limit  the  growth  of  numerical  instabilities  without  severely 
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compromising  the  accuracy  of  the  exospheric  transport  model.  The  adjustment  of  this  truncation 
parameter,  as  well  as  the  characteristic  exobase  radius,  R,  are  left  as  tasks  for  future  work. 


3.2.4  Software  Development 

Table  1  gives  an  overview  of  TIE-GCM-He  version  evolution  and  capabilities 


Table  2:  Software  development  history. 

Version/Date 

VI  .0 

Description  of  Commit 

First  version  to  run  with  4  major  species  02,  O,  He,  N2,  still  need  to  make  several  changes  throughout 

12/06/2012 

the  model.  Any  instance  of  N2  or  mbar  being  calculated,  need  to  be  changed  from  1-o2-o1  to  1-o2-o1-he 
and  from  1/(o2/mo2+o1/mo1+n2/mn2)  to  1/(o2/mo2+o1/mo1+he/mhe+n2/mn2). 

VI. 1 

12/11/2012 

Updated  addiag.F  so  that  Helium  is  taken  into  account  when  mbar,  xnmbar,  xnmbari,  xnmbarm,  DEN, 
Z,  ZG  and  others  are  calculated.  Changed  small  value  threshold  from  1e-6  to  1e-9  for  Helium  only. 
Changed  the  section  that  insures  non-negative  n2  values:  previously  the  denominator  differed  when 
normalizing  ol  and  o2  (this  seems  to  be  a  bug  in  the  original  code). 

VI  .2 

Added  thermal  diffusion  to  comp.F  for  Helium. 

12/12/2012 

VI  .3 

12/13/2012 

Updated  cpktkm.F  to  take  Helium  into  account  when  calculating  specific  heat,  molecular  diffusion  of  heat 
(conduction)  and  molecular  viscosity. 

VI  .4 

06/03/2013 

Modified  the  calculation  of  local  variable  xnmbar  in  comp_o2o.F  to  account  for  Helium;  the  variable  is 
used  to  calculate  prod/loss  rates  of  o2,  ol,  and  he.  Floor(sflux)  in  qrj.F  is  changed  from  80%  to  10% 
of  sfmin,  this  only  comes  into  play  during  very  low  solar  flux;  per  Stan  Solomons  instruction  (see  email 
dated  17-Jan-2013).  Made  co2u  variable  within  subroutine  newton  time-varying  to  capture  the  change  in 
C02  over  time  in  order  to  make  a  run  that  is  consistent  with  Solomon  et  al.  201 1 . 

VI  .5 

07/01/2013 

Specified  a  diffusive  flux  upper  boundary  condition  in  the  major  species  equation  describing  the  lateral 
exospheric  transport  of  helium  occurring  above  the  model  upper  boundary.  This  required  substantial 
modification  of  comp.F,  as  well  as  a  very  slight  modification  of  cons.F,  cons_variableC02.F,  and  dynam¬ 
ics.  F. 

VI  .6 

08/19/2013 

This  is  the  last  internal  AFRL  version  separate  from  NCAR/HAO  resulting  from  meeting  with  NCAR 
(S.  Solomon,  W.  Wang,  and  B.  Foster).  Identified  several  instances  of  (1-o2-o1),  and  implemented  a 
new  variable  n2=(1-o2-o1-he)  to  fix  all  occurrences  throughout  model.  Reestablished  MPI  functionality. 
Inducted  model  into  the  NCAR  software  repository  as  a  branch  of  TIE-GCM  to  facilitate  future  model 
development. 

NCAR/HAOs  subversion  revision  numbers  are  adopted  hereafter 

r966 

08/19/2013 

tiegcm/branches/tiegcm_he: 

-  Merge  this  branch  from  its  origin  at  tag  tiegcml  .94.2  with  the  current  trunk  (r950). 

r967 

08/20/2013 

tiegcm/branches/tiegcm_he: 

-  Added  N2  as  a  new  3d  field  in  fields. F  (defined  in  addiag.F). 

-  Replaced  code  using  (1  .-o2-o1  -he)  with  the  new  N2  field  (addiag.F,  comp.F,  comp_o2o.F,  and  cpktkm.F). 
N2  is  passed  to  these  routines  via  the  calls  in  advance  and  dynamics. 

-  At  this  point,  the  code  produces  identical  results  as  MPI  runs  from  the  previous  commit,  and  is  therefore 
still  very  close  to  Eric’s  original  non-mpi  code. 

-  Next,  i  will  replace  code  using  (1  ,-o2-o1 )  with  n2  (which  is  now  1  .-o2-o1  -he),  so  after  that  is  committed, 
results  will  be  different  than  Eric’s  original  code.  This  occurs  in  several  routines  (chapman,  comp_n4s, 
comp no,  dt,  duv,  elden,  lamdas,  minor,  etc) 

r968 

08/21/2013 

tiegcm/branches/tiegcm_he: 

-  Replaced  all  occurrences  of  1  ,-o2-o1  with  the  new  n2  3d  array,  which  is  now  1-o2-o1-he.  This  means 
the  results  will  differ  from  the  previous  commit,  and  Eric’s  original  non-mpi  code.  Results  are  presumably 
more  accurate  because  now  helium  is  included  in  N2. 

-  Added  new  N2  diagnostic  in  diags.F.  mkdiag_N2  is  called  from  comp.F  to  save  updated  N2  to  secondary 
histories.  Now  users  can  add  ’N2’  to  namelist  SECFLDS  without  having  to  call  addfld. 

-  Note  helium  was  not  included  in  the  calculation  of  the  bottom  interface  level  of  n2  (n2i(1 ,..))  in  qinite.F 
and  qrj.F.  For  now,  1  am  using  n2i(lev0,..)  =  1.-o2i(lev0,..)-o1i(lev0,..)  because  helium  is  negligible  at  this 
bottom  boundary  level.  Search  on  ”btf  tiegcm.he”  to  find  this  in  the  source  code. 

r973 

08/28/2013 

tiegcm/branches/tiegcm_he: 

-  The  problem  described  in  the  previous  revision  (r969)  was  a  result  of  a  misplaced  parenthesis  in  rrk. 
This  commit  fixes  that  problem,  and  also  adds  He/0+  collision  frequency  provided  by  Wenbin. 
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rl  01 2 

02/03/2014 

tiegcm/branches/tiegcm_he: 

-  Update  this  branch  from  the  trunk,  including  all  changes  to  the  trunk  since  this  branch  was  started 
(r964).  This  was  done  w/  the  following  merge  command,  executed  from  an  up-to-date  working  copy  of 
the  branch:  svn  merge  -r964:HEAD  $SVN/tiegcm/trunk.  There  were  no  conflicts.  This  is  in  preparation 
for  merging  the  helium  branch  to  the  trunk. 

rl  01 3 

02/03/2014 

tiegcm/trunk: 

-  Merge  the  tiegcm  branch  tiegcm.he  to  the  trunk.  This  was  done  with  the  following  command,  executed 
from  an  up-to-date  working  copy  of  the  trunk: 

svn  merge  $SVN/tiegcm/trunk  $SVN/tiegcm/branches/tiegcm_he 

There  were  no  conflicts  with  this  merge,  which  follows  a  recent  merge  to  update  the  helium  branch  from 
the  trunk  (see  -rl  012). 

-  This  adds  helium  ’HE’  as  a  prognostic  in  the  model,  so  it  is  included  in  primary  histories,  along  with 
helium  at  the  previous  timestep  ’HE_NM’.  The  code  to  solve  for  helium  was  provided  by  Eric  Sutton, 
and  parallelized  by  myself.  See  especially  new  source  files  laplacian.F,  and  sh.coef.F.  All  calculations  of 
neutral  density,  and  N2  were  modified  appropriately.  N2  was  added  as  a  3d  diagnostic. 

-  Many  of  Eric’s  comments  (’’Added  by  EKS”,  "Modified  by  EKS”,  etc)  were  removed  for  this  commit  to 
the  trunk  (all  of  them  were  retained  in  the  tiegcm.he  branch). 

-  For  more  details  regarding  the  addition  of  helium  in  the  model,  please  see  the  svn  commit  log  for 
$SVN/tiegcm/branches/tiegcm he. 

rl  047 

05/08/2014 

tiegcm/trunk: 

-  Helium  at  2.5-degree  resolution  is  now  working,  thanks  to  Eric  Sutton  for  providing  laplacian  coeffi¬ 
cients  for  this  resolution.  See  new  source  files  he_coefO_dres.F  (Ipmn),  and  he_coef1_dres.F  (zmn).  Also 
changed  the  name  of  coefficients  file  for  5-deg  resolution  from  sh_coef.F  to  he_coefs_sres.F.  These  three 
files  are  modules  with  data  statements  dependent  on  resolution  as  determined  by  cpp  #if  conditionals. 
They  are  use-associated  by  laplacian.F. 

-  Set  calc.helium  namelist  flag  to  1  for  both  resolution  default  namelist  input  files  in  the  scripts  directory 
(see  r1042).  Removed  trap  to  shutdown  if  calc helium=1  at  2.5-degree  resolution. 

rl  061 

06/23/2014 

tiegcm/trunk: 

-  Set  HE  lower  boundary  (mmr)  pshelb  =  0.1154E-5  (see  Ibc.F)  (in  timegcm,  pshelb=7.18e-7).  Pshelb  is 
used  in  subs  bndef  and  bndcmp  (Ibc.F)  to  set  the  lower  boundary  (gg,  fb,  etc).  If  a  history  is  read  that 
does  NOT  have  helium,  then  the  3d  (global  and  vertical)  HE  and  HE_NM  are  initialized  to  pshelb  (see 
rdsource.F).  This  helps  "prime”  helium  on  its  approach  to  steady  state  without  a  starting  field  from  the 
history.  Continuation  runs  of  course  read  helium  from  the  given  history  from  which  to  continue. 

-  A  5-deg  run  of  tiegcm  from  March  equinox  solar  minimum  conditions,  showed  that  it  can  take  80  days 
for  helium  to  reach  approximate  steady  state,  starting  from  global  pshelb. 

1  will  be  making  a  similar  commit  in  the  timegcm  trunk. 

rl  060 

06/16/2014 

tiegcm/trunk: 

-  Pass  helium  to  minor  through  comprno  and  comp_n4s,  so  it  can  be  included  in  xmbari,  xmbar_kh,  etc. 
See  minor.F. 

-  To  be  completely  correct  and  consistent,  we  will  also  have  to  add  pshe,  dhedz  (analogous  to 
pso2,do2dz,  etc),  and  make  alfa  3x3  instead  of  2x2  so  minor  species  diffusion  through  the  major  species 
will  include  helium  (see  ex(k,i)).  See  ”btf  6/16/14”  comment  in  minor.F.  1  am  working  with  Wenbin  on  this. 

rl  074 

08/14/2014 

tiegcm/trunk: 

-  Add  argon  as  a  minor  species.  This  is  similar  to  the  r1069  commit  which  added  argon  to  the  timegcm 
trunk.  New  source  file  comp.ar.F. 

-  As  in  timegcm,  global  mean  column  values  from  glbmean  model  are  used  to  initialize  argon  when  it  is 
not  available  on  the  source  history. 

-  Note  that  tiegcm  does  not  calculate  horizontal  diffusion  for  minor  species,  like  timegcm  does. 

-  Test  runs  have  been  made  at  both  5.0-deg  and  2. 5-deg  resolutions,  but  further  testing  and  validation  of 
argon  should  be  done  with  this  code  before  the  next  major  release  (v2.0). 
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rl  1 1 9 

01/21/2015 

tiegcm/trunk  (helium): 

-  Modifications  to  helium  calculation  recommended  by  Eric  Sutton, during  his  visit  to  hao  this  week  (see 
also  email  correspondence  in  Nov,  2014): 

-  Deleted  data  array  Ipmn  which  contained  (m''2/sin(theta)-n*(n+1))*pmn.  The  m~2/sin(theta)  term  is  not 
needed.  Therefore,  1  replaced  Ipmn  with  pmn  (the  normalized  associated  legendre  polynomials)  and 
added  the  latitude-independent  term  -n*(n+1)  to  the  synthesis  loop  within  laplacian.F.  (src/coefs_sres.F 
and  src/coefO_dres.F)  (makes  this  commit  a  large  change  set) 

-  Added  compiler  (preprocessor)  #IF  statement  and  parameter  ’truncdeg’  to  define  the  degree  of  trunca¬ 
tion  for  the  spectral  approximation,  which  depends  on  grid  resolution  and  upper  boundary  pressure  level. 
Set  truncdeg=4  for  default  (5-deg)  resolution,  and  truncdeg=8  for  double  (2.5-deg)  resolution. 

-  Minor  updates  to  comments  in  laplacian.F,  comp.F,  he_coefs_sres.F,  and  he_coef0_dres.F. 

-  The  ’Fit’  and  ’Synthesis’  sections  of  laplacian.F  each  now  contain  commented  loops  (just  after  the  main 
loops)  that  are  helpful  when  debugging  various  levels  of  truncation  and  testing  for  stability. 

-  Removed  unneeded  z  argument  in  comp.F  (and  removed  it  from  the  call  to  comp  in  dynamics. F). 

-  Left  addfld  calls  for  FLX HE  and  FLX ARG  uncommented  in  comp.F  (after  call  to  laplacian). 

rl  125 

01/28/2015 

tiegcm/trunk: 

-  The  following  comments  are  from  Eric  Sutton  (see  correspondence  1/22/15): 

-  A  temperature  dependent  variation  was  applied  to  the  coefficients  describing  diffusion  of  helium  with 
the  other  major  species  (comp.F).  In  previous  versions,  the  temperature  and  pressure  dependence  of 
diffusion  coefficients  were  normalized  by  assuming  a  variation  of  the  form  (p_00/p)*(T/T_00),'1 .75.  In  the 
case  of  helium,  the  exponent  1 .75  is  not  realistic,  requiring  a  correction  factor  of  (T/T_00)"(1 .75-s),  where 
s  is  the  appropriate  exponent,  to  be  multiplied  to  the  normalized  diffusion  coefficients  (phi’s)  describing 
He-02,  He-O,  and  He-N2  diffusion.  The  correction  factor  is  generally  on  the  order  of  5%. 

3.3  Results  and  Discussion 

In  this  section,  we  present  the  salient  features  of  the  new  model.  While  many  simulations 
were  necessary  in  order  to  distill  our  description  of  these  features  with  respect  to  season,  local 
time,  latitude,  external  forcing  parameters,  and  boundary  conditions,  only  a  small  subset  of  sim¬ 
ulations  are  presented.  These  were  created  using  the  model  settings  and  inputs  described  in 
Section  ^T27R  and  are  specific  to  the  prevailing  solar  and  geophysical  conditions  of  2008.  Sec¬ 
tion^]  includes  additional  plots  to  aid  in  visualization,  specifically  regarding  the  sensitivity  of  the 
helium  distribution  to  external  forcing  and  boundary  conditions. 

Figure  3J] shows  helium  densities  at  250  km  altitude  simulated  by  the  TIE-GCM  during  each 
of  the  four  seasons  of  2008.  The  winter  helium  bulge  phenomenon  is  clearly  present  at  both 
solstices.  During  the  equinoxes,  the  helium  bulge  undergoes  a  migration  from  the  spring  hemi¬ 
sphere  to  the  fall.  Along  the  way,  helium  levels  are  briefly  enhanced  at  low  latitudes  with  a 
strong  preference  for  early  morning  local  times,  with  the  full  transistion  taking  approximately  1-2 
months.  At  winter  solstice,  a  similar  preference  for  early  morning  is  tempered  by  an  aversion 
to  the  auroral  zones,  where  pockets  of  divergence  and  upwelling  lead  to  localized  helium  de¬ 
pletions.  This  balance  manifests  as  a  diurnal  modulation  of  the  winter  helium  bulge  in  latitude 
and  local  time.  Symptoms  of  this  behavior  can  be  seen  in  the  the  upper  right  panel  of  Fig¬ 
ured  where  the  southern  hemisphere  winter  peak  occurs  around  16:00  LT.  For  reference,  the 
geomagnetic  poles  are  located  at  82.4°N/18:30  LT  and  74.5°S/8:20  LT  in  these  plots. 

The  high-latitude  helium  distribution  is  further  complicated  by  short-scale  variations  in  geo¬ 
magnetic  heating.  In  general,  helium  densities  tend  to  increase  at  low  latitudes  during  periods  of 
geomagnetic  activity.  The  opposite  is  true  in  the  polar  region  during  solstice,  as  the  high-latitude 
upwelling  and  divergence  resulting  from  geomagnetic  activity  tend  to  lift  heavy  constituents  while 
dispersing  helium  over  a  larger  horizontal  expanse. 

The  distribution  of  helium  is  highly  sensitive  to  geomagnetic  activity,  the  effects  of  which  can 


25 

Approved  for  public  release;  distribution  is  unlimited. 


cm  3  xIO6  cm  3  x  107 


Figure  11:  Helium  concentration  (TIE-GCM).  Global  distribution  of  helium  number  densities 
at  250  km  altitude  during  each  season  for  solar  minimum  conditions  (2008),  as  calculated  by 
TIE-GCM.  Equinox  plots  (left)  share  a  common  color  scale,  as  do  solstice  plots  (right). 


be  seen  in  the  contrasting  equinoctial  helium  distributions  of  Figure  HI  The  March  equinox  con¬ 
sists  of  enhanced  low  and  middle  latitude  helium  densities  accompanied  by  depletions  closer  to 
the  poles,  all  associated  with  a  slight  elevation  in  the  level  of  geomagnetic  activity  over  the  pre¬ 
vious  3-hour  period  (KP=2.0)  relative  to  the  September  equinox  (KP=0.3).  The  same  argument 
can  be  applied  to  the  solstice  plots  of  Figure  [Til  wherein  the  slightly  disturbed  (KP=2.0)  June 
solstice  helium  distribution  is  shifted  away  from  the  winter  pole  in  comparison  to  the  undisturbed 
(KP=0.0)  December  solstice.  The  helium  distribution  is  most  certainly  influenced  by  the  time  his¬ 
tory  of  geomagnetic  activity  over  the  previous  ~24  hours  or  more.  As  such,  an  index  describing 
the  level  of  geomagnetic  activity  over  a  3-hour  interval  may  not  generally  be  a  reliable  indicator. 
However,  in  all  four  of  the  cases  presented  the  3-hour  KP  index  is  fairly  representative  of  the 
levels  of  geomagnetic  activity  during  the  previous  24-hour  period.  The  solstice  comparison  is 
less  straightforward  than  for  equinox  due  to  several  additional  complications.  One  such  compli¬ 
cation  is  that  the  location  of  maximum  helium  concentration  is  more  sensitive  during  solstice  to 
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the  location  of  the  geomagnetic  poles.  The  solstice  comparisons  also  suffer  from  slightly  differ¬ 
ing  amounts  of  solar  flux.  The  Supporting  Information  provides  additional  figures  emphasizing 
the  sensitivity  of  the  helium  distribution  throughout  the  year  to  variations  in  geomagnetic  activity, 
solar  flux,  and  forcing  of  the  lower  boundary  by  migrating  tides. 

As  a  basis  for  comparison,  Figure  32]  shows  helium  densities  at  250  km  altitude  as  calcu¬ 
lated  by  the  MSIS  model.  Many  of  the  salient  features  are  qualitatively  similar  to  those  of  the 
TIE-GCM,  with  respect  to  seasonal,  latitudinal,  and  local  time  characteristics.  MSIS  helium  dis¬ 
tributions  clearly  exhibit  the  same  strong  preference  for  the  winter  polar  regions  during  solstice, 
and  for  the  low-latitude,  early  local  time  sectors  during  equinox.  Likewise,  a  similar  sensitivity  to 
geomagnetic  effects  is  evident  within  MSIS.  Notice,  however,  that  the  color  scales  differ  between 
Figures  33] and  3H in  order  to  show  behavior  over  the  full  range  of  each  model.  At  250  km,  the 
TIE-GCM  typically  underestimates  the  magnitude  of  the  MSIS  helium  bulge  by  approximately 
20%  during  solstice,  while  overestimating  it  by  5%  during  equinox.  This  agreement  is  reason¬ 
able,  considering  that  no  adjustments  have  been  made  to  the  TIE-GCM  in  an  effort  to  improve 
model  agreement. 

Likewise,  the  MSIS  model  estimated  and  applied  correction  factors  for  the  underlying  mass 
spectrometer  data  [13],  which  could  further  limit  the  absolute  accuracy  of  such  model  compar¬ 
isons.  In  certain  cases,  there  are  discrepancies  in  the  location  and  shape  of  the  helium  bulge 
between  models.  For  instance,  the  location  of  maximum  helium  concentration  during  the  June 
solstice  is  out  of  phase  by  about  8  hours  in  local  time  between  the  two  models.  While  the  MSIS 
helium  distribution  is  prescribed,  to  a  certain  extent,  by  a  trade-off  between  the  data  sparsity  of 
its  underlying  historical  data  set  and  the  complexity  of  its  basis  functions,  further  investigation  is 
needed  before  attributing  any  discrepancies  to  the  shortcomings  of  either  model. 

Figure  33]  shows  the  magnitude  of  the  helium  bulge  ratio  as  a  function  of  height,  during  solar 
minimum  solstice  conditions.  These  profiles  were  constructed  by  taking  the  ratio  of  maximum- 
to-minimum  helium  number  densities  along  each  model  meridian  to  roughly  approximate  the 
method  of  calculation  used  in  previous  studies.  The  ratio  at  each  height  was  then  averaged 
both  zonally  and  over  the  course  of  a  day;  note  that  no  attempt  was  made  to  specify  the  local 
time  sampling  of  a  particular  polar-orbiting  satellite.  The  vertical  profiles  exhibit  a  quick  increase 
from  the  lower  boundary,  giving  way  to  a  maximum  around  175  km,  then  decaying  slowly  with 
altitude  to  the  upper  boundary.  This  behavior  can  be  explained  by  the  transition  from  a  region 
below  the  peak  which  is  dominated  by  collisions,  to  a  region  above  the  peak  where  diffusive 
equilibrium  is  well  established.  Below  the  height  of  maximum  bulge  ratio,  the  summer-to-winter 
bulk  circulation  pattern  leads  to  the  accumulation  of  helium  in  the  winter  hemisphere.  Above 
this  height,  however,  vertical  profiles  begin  to  approximate  diffusive  equilibrium,  causing  helium 
densities  in  the  winter  hemisphere  to  decrease  with  height  at  a  slightly  faster  exponential  rate 
than  those  in  the  warmer  summer  hemisphere. 

The  significant  difference  between  June  and  December  is  due  to  a  combination  of  lower  solar 
flux  and  geomagnetic  activity  during  the  December  solstice.  Smaller  contributions  to  this  differ¬ 
ence  may  arise  from  seasonal  variations  such  as  in  the  eddy  diffusivity.  Error  bars  in  Figure  331 
show  the  standard  deviation  of  the  helium  ratio  over  the  course  of  a  day,  giving  an  indication 
of  the  sensitivity  to  diurnal  variations  as  well  as  small  variations  in  geomagnetic  activity.  Below 
150  km,  smaller  standard  deviations  are  seen,  indicating  that  variations  in  the  lower  part  of  the 
profile  take  place  on  longer  timescales.  Presumably,  the  lower  portion  of  the  profile  is  more  sen¬ 
sitive  to  season  and  solar  flux  than  to  short-scale  geomagnetic  activity.  Approaching  altitudes  as 
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Figure  12:  Helium  concentration  (MSIS).  Global  distribution  of  helium  number  densities  at 
250  km  altitude  during  each  season  for  solar  minimum  conditions  (2008),  as  calculated  by  MSIS. 
Equinox  plots  (left)  share  a  common  color  scale,  as  do  solstice  plots  (right);  these  are  distinct 
from  the  color  scales  of  Figure  11'. 


low  as  100  km,  the  two  profiles  begin  to  converge,  suggesting  a  muted  response  to  geomagnetic 
activity  as  well  as  to  seasonal  variations. 

The  addition  of  helium  to  the  TIE-GCM  has  several  feedback  effects  on  the  global  structure  of 
the  model.  Most  of  these  are  related  to  the  change  in  the  mean  mass,  which  can  become  quite 
small  and  even  approach  4  amu  near  the  top  of  the  model.  On  levels  of  constant  pressure,  such 
a  decrease  in  the  mean  mass  corresponds  directly  to  a  decrease  in  mass  density.  At  a  fixed 
height,  however,  this  behavior  is  accompanied  by  the  expansion  of  the  atmosphere  according  to 
the  ideal  gas  law,  causing  levels  of  constant  pressure  to  move  upward.  With  increasing  altitude, 
the  expansion  effect  begins  to  dominate  the  mean-mass  effect  such  that  the  decay  in  mass 
density  with  height  becomes  much  more  gradual  when  helium  is  considered.  Figure  T41  shows 
the  induced  increase  in  mass  density  at  a  fixed  altitude  of  415  km.  While  the  inclusion  of  helium 
causes  the  model’s  upper  boundary  to  expand  considerably  higher  than  415  km,  we  chose  this 
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Figure  13:  Winter  helium  bulge  ratio.  Vertical  profile  of  the  winter-to-summer  helium  bulge 
ratio  during  solar  minimum  June  (black,  June  21,  2008)  and  December  (grey,  Dec.  21,  2008) 
solstice  conditions.  The  profiles  represent  the  daily  average  of  the  ratio  of  maximum-to-minimum 
helium  number  densities  taken  along  each  meridian,  roughly  approximating  the  sampling  of  a 
polar  orbiting  satellite  (see  text  for  a  detailed  explanation).  Error  bars  indicate  the  standard 
deviation  of  values  over  the  course  of  a  day. 


height  for  our  comparison  because  it  was  the  highest  altitude  that  remained  within  the  vertical 
domain  of  the  original  TIE-GCM  simulations  during  each  of  the  four  time  periods  shown. 

The  increase  in  mass  density  is  most  noticeable  during  solstice,  where  differences  of  20-25% 
can  be  seen.  Both  equinox  and  solstice  mass  density  increases  are  largest  under  quiet  geomag¬ 
netic  conditions.  While  somewhat  modest,  these  percent  differences  increase  with  height  at  an 
approximate  rate  of  1  %  per  kilometer  near  the  upper  boundary  of  the  TIE-GCM  in  regions  of  large 
helium  densities.  If  the  composition  of  the  TIE-GCM  is  extended  vertically  into  the  exosphere 
under  the  first-order  approximation  of  diffusive  equilibrium,  the  effects  of  helium  soon  become 
the  dominant  factor  in  neutral  mass  density  variations.  Under  solar  minimum  conditions,  an  ex¬ 
tension  of  both  models  to  500  km  results  in  differences  on  the  order  of  50%  during  equinox  and 
100-200%  during  solstice.  At  600  km,  the  solstice  differences  exceed  of  an  order  of  magnitude. 

Contours  in  Figure25]show  differences  in  the  height  of  a  log-pressure  level  near  the  top  of  the 
model  induced  by  the  inclusion  of  helium.  Near  the  winter  pole  where  these  height  differences 
maximize,  the  atmosphere  is  uplifted  by  some  50-60  km  when  compared  to  an  atmosphere  sim¬ 
ulated  without  helium.  This  modification  further  couples  to  the  horizontal  momentum  equations 
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Figure  14:  Helium’s  impact  on  mass  density.  Percent  increase  in  the  total  mass  density  at  a 
fixed  altitude  of  415  km  resulting  from  the  inclusion  of  helium  in  TIE-GCM  during  each  season 
for  solar  minimum  conditions  (2008).  Equinox  plots  (left)  share  a  common  color  scale,  as  do 
solstice  plots  (right). 


[see  0],  increasing  horizontal  gradients  in  the  geopotential  and  resulting  in  a  difference  wind  pat¬ 
tern  that  flows  away  from  the  winter  helium  bulge,  as  depicted  by  the  vector  arrows  of  Figure  451 
This  effect  generally  becomes  noticeable  in  the  upper  thermosphere,  above  300-400  km,  where 
differences  as  high  as  15-20  m/s  can  be  attained. 

3.4  Conclusions 

This  chapter  establishes  methods  for  tracking  helium  abundance  self-consistently  through¬ 
out  the  thermosphere.  The  resulting  model  simulations  qualitatively  recreate  the  expected  sea¬ 
sonal/latitudinal  behavior  while  also  showing  reasonable  quantitative  agreement  with  MSIS.  More¬ 
over,  the  model  provides  winter-to-summer  helium  ratios  that  generally  agree  with  solar  minimum 
observations  from  AE-C  [[2J.  A  more  rigorous  one-to-one  comparison  between  this  new  model 
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Figure  15:  Helium’s  impact  on  geopotential  height  and  horizontal  winds.  Difference  in  the 
geopotential  height  (color  contours)  and  the  horizontal  wind  field  (vectors)  on  a  level  of  constant 
pressure  near  the  upper  model  boundary  (^=+6.75)  resulting  from  the  inclusion  of  helium  in 
TIE-GCM  during  each  season  for  solar  minimum  conditions  (2008).  Equinox  plots  (left)  share 
common  color  and  vector  scales,  as  do  solstice  plots  (right). 
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and  legacy  mass  spectrometer  measurements  is  merited;  however,  this  task  is  left  for  future 
work. 

Perhaps  the  most  direct  application  for  this  new  model  is  related  to  the  increased  realism  of 
the  neutral  mass  density  vertical  profile,  and  thus  the  improvement  in  model  performance  with 
respect  to  satellite  drag  observations  in  the  upper  thermosphere.  At  a  constant  height  within 
the  model  domain,  we  have  shown  that  including  helium  in  the  TIE-GCM  causes  differences 
in  neutral  mass  density  on  the  order  of  20-30%  during  solar  minimum.  The  most  noticeable 
differences  occur  near  the  upper  model  boundary  during  solstice  in  the  winter  hemisphere. 

Furthermore,  helium  concentration  in  the  exosphere  is  highly  sensitive  to  the  dynamics  of 
the  thermosphere.  An  appropriate  exospheric  model  could  use  the  TIE-GCM’s  upper  boundary 
to  specify  a  realistic  exobase.  Using  profiles  approximated  by  diffusive  equilibrium  above  the 
TIE-GCM’s  upper  boundary,  we  demonstrated  that  helium  can  account  for  order-of-magnitude 
differences  in  neutral  density  near  600  km  and  above.  These  differences,  structured  in  latitude 
and  local  time,  are  strongly  modulated  by  season  and  geomagnetic  activity,  lending  significant 
variability  to  the  upper  thermosphere  and  exosphere.  This  seasonal,  latitudinal,  and  local  time 
helium  behavior  can  be  used  to  inform  the  structure  of  semi-empirical  model  basis  functions  [e.g. 
71].  At  a  minimum,  inferring  the  amplitude  of  such  basis  functions  would  require  sufficient  cover¬ 
age  of  high-altitude  satellite  drag  measurements,  but  would  be  better  served  by  a  contemporary 
set  of  mass  spectrometer  measurements. 

The  value  of  helium  as  a  tracer  of  thermospheric  dynamics  has  been  known  for  some  time 
[see  (52[].  In  addition  to  its  ability  to  diagnose  the  interplay  of  circulation  and  diffusion  in  the 
thermosphere,  our  new  model  will  enable  future  studies  attempting  to  exploit  the  sensitivity  of 
the  helium  distribution  to  otherwise  unobservable  system  dynamics  and  inputs.  We  anticipate 
that  employing  helium  as  a  diagnostic  tracer — e.g.  in  order  to  specify  or  constrain  high-latitude 
energy  inputs,  solar-driven  circulation  pattern  strength,  and/or  sub-grid  scale  model  dynamics — 
will  be  beneficial  in  refining  model  performance  for  scientific  endeavors  as  well  as  operational 
applications. 


3.5  Additional  Supporting  Information 

The  supporting  information  contains  several  controlled  runs  illustrating  the  sensitivity  of  the 
helium  distribution  to  geomagnetic  activity,  solar  flux,  and  migrating  tides  imposed  at  the  lower 
boundary. 
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Figure  16:  TIE-GCM  control  run:  constant  M10.7  and  KP.  Global  distribution  of  helium  number 
densities  at  250  km  altitude  during  each  season.  Solar  and  geomagnetic  drivers  were  held 
constant  throughout  the  year-long  simulation  at  the  median  values  of  2008:  M10.7  =  62.0  and 
KP  =  1.3.  Equinox  plots  (left)  share  a  common  color  scale,  as  do  solstice  plots  (right). 


33 

Approved  for  public  release;  distribution  is  unlimited. 


0  4  8  12  16  20  24  0  4  8  12  16  20  24 

Local  Time  (hours)  Local  Time  (hours) 


cm  3  x10G  cm  3  x  107 


Figure  17:  TIE-GCM  control  run:  elevated  KP.  Global  distribution  of  helium  number  densities 
at  250  km  altitude  during  each  season.  Solar  and  geomagnetic  drivers  were  held  constant  as 
in  Figure  3H  but  with  an  elevated  geomagnetic  index  of  KP  =  3.0.  Equinox  plots  (left)  share  a 
common  color  scale,  as  do  solstice  plots  (right). 
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Figure  18:  TIE-GCM  control  run:  elevated  M10.7.  Global  distribution  of  helium  number  densities 
at  250  km  altitude  during  each  season.  Solar  and  geomagnetic  drivers  were  held  constant  as 
in  Figure  3H  but  with  an  elevated  solar  flux  proxy  of  M10.7  =  100.0.  Equinox  plots  (left)  share  a 
common  color  scale,  as  do  solstice  plots  (right). 


35 

Approved  for  public  release;  distribution  is  unlimited. 


Local  Time  (hours)  Local  Time  (hours) 


Figure  19:  TIE-GCM  control  run:  Tidal  Forcing  Removed.  Global  distribution  of  helium  num¬ 
ber  densities  at  250  km  altitude  during  each  season.  Solar  and  geomagnetic  drivers  were  held 
constant  as  in  Figure  HU  but  with  the  lower  boundary  GSWM  forcing  from  migrating  tides  re¬ 
moved.  Equinox  plots  (left)  share  a  common  color  scale,  as  do  solstice  plots  (right). 
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A.1  Final  Laboratory  Task  Report 

LRIR  #:  13RV09COR 

Reporting  Period:  Fiscal  Year  2013-2015 

Laboratory  Task  Manager:  Eric  Sutton 

Commercial  Phone:  (505)  846-7846  DSN:  246-7846 

Mailing  Address: 

RVBXI  /  Bldg.  570  /  Rm.  2363 
3550  Aberdeen  Ave  SE 
Kirtland  AFB,  NM  87117 

E-Mail  Address:  eric.sutton.5@us.af.mil 

AFOSR  Program  Manager:  Kent  Miller 


Table  A-1:  Funding  Summary  by  Cost  Category 


Year 

In-House 

Capital  Equip. 
(>$  5K  each) 

Non-Capital  Equip. 
(<$  5K  each) 

Travel 

On-Site 

Contractor 

Total 

FY13 

138K 

— 

— 

3K 

9K 

150K 

FY14 

142K 

— 

3K 

5K 

— 

150K 

FY15 

144K 

— 

— 

6K 

— 

150K 

Table  A-2:  Personnel 


Name 

Degree 

Discipline 

Involvement 

Air  Force  Employees: 

Eric  Sutton 

PhD 

Physicist 

0.38  Man-Yr/Yr 

Samuel  Cable 

PhD 

Physicist 

0.10  Man-Yr/Yr 

A.2  Invited  Talks  and  Lectures 

Sutton,  E.  K.  (2015),  Dynamical  influences  on  thermospheric  composition  (invited),  AGU  Fall 
Meeting,  Dec.  2015  (tentatively  scheduled). 

Sutton,  E.  K.,  and  F.  Schmidt  (201 5),  First  Principle  Modeling  of  the  Upper  Thermosphere/Lower 
Exosphere  (invited),  CEDAR  Workshop,  Seattle,  WA,  21-25  June,  2015. 
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Sutton,  E.  K.  (2015),  Satellite  Accelerometer  Observations  During  the  Storms  and  Substorms 
Without  Borders  Campaigns  (invited),  CEDAR  Workshop,  Seattle,  WA,  21-25  June,  2015. 

Sutton,  E.  K.  (2015),  The  Upper  Atmosphere  and  its  Influence  on  Satellite  Motion  (invited  Lec¬ 
ture),  University  of  New  Mexico,  ECE  graduate  course  on  Space  Weather,  Albuquerque,  NM, 
27  April,  2015. 

Sutton,  E.  K.  (2013),  The  Role  of  Helium  in  the  Thermosphere  During  Recent  Solar  Minima, 
Invited  Guest  Lecturer,  HAO  Colloquium  Series,  NCAR/HAO,  Boulder,  CO,  8  August  2013. 

Sutton,  E.  K.  (2011),  Upper  Atmosphere  Modeling  at  AFRL,  Invited  Speaker,  1st  IMPACT  Team 
Meeting,  LANL,  Santa  Fe,  NM,  19  October,  2011. 


A.3  Scientific  Visits 

Sutton,  E.  K.  (2015),  NCAR/HAO  Visiting  Scientist  Program,  Boulder,  CO,  19-21  January  2015. 
Sutton,  E.  K.  (2013),  NCAR/HAO  Visiting  Scientist  Program,  Boulder,  CO,  5-9  August  2013. 


A.4  Publications 

A.4.1  First-Authored 

Sutton,  E.  K.,  J.  P.  Thayer,  W.  Wang,  S.  C.  Solomon,  X.  Liu,  and  B.  T.  Foster  (2015),  A  self- 
consistent  model  of  helium  in  the  thermosphere,  Journal  of  Geophysical  Research:  Space 
Physics,  pp.  n/a-n/a,  doi:1 0.1 002/201 5JA021 223. 

Sutton,  E.  K.,  S.  B.  Cable,  C.  S.  Lin,  L.  Qian,  and  D.  R.  Weimer  (2012),  Thermospheric  ba¬ 
sis  functions  for  improved  dynamic  calibration  of  semi-empirical  models,  Space  Weather,  10, 
SI 0001 ,  doi:1 0.1 029/201 2SW000827. 

A.4.2  Co-Authored 

Huang,  C.  Y.,  Y.  Huang,  Y.-J.  Su,  E.  K.  Sutton,  E.  K.,  M.  Hairston,  and  W.  R.  Coley  (2015,  ac¬ 
cepted),  Ionosphere-Thermosphere  Response  to  Solar  Wind  Forcing  During  Magnetic  Storms, 
accepted  for  publication  in  the  J.  Space  Wea.  Space  dim. 

Yamazaki,  Y.,  M.  J.  Kosch,  and  E.  K.  Sutton  (2015),  North-south  asymmetry  of  the  high-latitude 
thermospheric  density:  IMF  BY  effect,  Geophys.  Res.  Lett.,  42,  225-232, 
doi:1 0.1 002/201 4GL062748. 

Liu,  X.,  W.  Wang,  J.  P.  Thayer,  A.  Burns,  E.  K.  Sutton,  S.  C.  Solomon,  L.  Qian,  and  G.  Lucas 
(2014),  The  winter  helium  bulge  revisited,  Geophys.  Res.  Lett.,  41, 6603-6609, 
doi:1 0.1 002/201 4GL061 471. 

Vickers,  H.,  M.  J.  Kosch,  E.  K.  Sutton,  L.  Bjoland,  Y.  Ogawa,  and  C.  La  Hoz  (2014),  A  solar  cycle 
of  upper  thermosphere  density  observations  from  the  EISCAT  Svalbard  Radar,  J.  Geophys. 
Res.  (Space  Physics),  119,  6833-6845,  doi:1 0.1 002/201 4JA01 9885. 

Liu,  X.,  J.  P.  Thayer,  A.  Burns,  W.  Wang,  and  E.  K.  Sutton  (2014),  Altitude  variations  in  the  ther¬ 
mosphere  mass  density  response  to  geomagnetic  activity  during  the  recent  solar  minimum,  J. 
Geophys.  Res.  (Space  Physics),  119,  2160-2177,  doi:1 0.1 002/201 3JA01 9453. 
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Huang,  C.  Y.,  Y.-J.  Su,  E.  K.  Sutton,  D.  R.  Weimer,  and  R.  L.  Davidson  (2014),  Energy  coupling 
during  the  August  2011  magnetic  storm,  J.  Geophys.  Res.  (Space  Physics),  119,  1219-1232, 
doi:1 0.1 002/201 3JA01 9297. 

Mehta,  P.  M.,  C.  A.  McLaughlin,  and  E.  K.  Sutton  (2013),  Drag  Coefficient  Modeling  for  GRACE 
using  Direct  Simulation  Monte  Carlo,  Adv.  Space  Res.,  52  (12),  doi:1 0.1 01 6/j.asr.201 3.08.033. 

Huang,  C.  Y.,  P.  A.  Roddy,  E.  K.  Sutton,  R.  Stoneback,  R.  F.  Pfaff,  L.  C.  Gentile,  S.  H.  Delay 
(2013),  Ion-neutral  coupling  during  deep  solar  minimum,  J.  Atmos.  Sol. -Terr.  Phys.,  103, 
doi:1 0.1 01 6/j.jastp.201 2.1 1.009. 

Vickers,  H.,  M.  J.  Kosch,  E.  K.  Sutton,  Y.  Ogawa,  and  C.  La  Hoz  (2013),  Thermospheric  atomic 
oxygen  density  estimates  using  the  EISCAT  Svalbard  Radar,  J.  Geophys.  Res.  (Space 
Physics),  118,  1319-1330,  doi :  1 0. 1 002/jgra.501 69 . 

Titov,  E.  V.,  J.  M.  Burt,  E.  Josyula,  E.  K.  Sutton,  and  C.  S.  Lin  (2013),  Uncertainties  in  Satellite 
Drag  Associated  with  Variations  in  Atmospheric  Conditions,  in  proc.  5th  AIAA  Atmospheric 
and  Space  Environments  Conference,  San  Diego,  CA,  24-27  June  2013,  AIAA  2013-2683, 
doi:10.251 4/6.201 3-2683. 

Mishin,  E.,  E.  Sutton,  G.  Milikh,  I.  Galkin,  C.  Roth,  and  M.  Frster  (2012),  F2-region  atmospheric 
gravity  waves  due  to  high-power  HF  heating  and  subauroral  polarization  streams,  Geophys. 
Res.  Lett.,  39,  L1 1 101 ,  doi:1 0.1 029/201 2GL052004. 

Lei,  J.,  T.  Matsuo,  X.  Dou,  E.  K.  Sutton,  and  X.  Luan  (2012),  Annual  and  semiannual  variations 
of  thermospheric  density:  EOF  analysis  of  CHAMP  and  GRACE  data,  J.  Geophys.  Res.,  117, 
A01 310,  doi:1 0.1 029/201 1 JA01 7324. 

Shim,  J.  S.,. . . ,  E.  K.  Sutton,  et  al.  (2012),  CEDAR  Electrodynamics  Thermosphere  Ionosphere 
(ETI)  Challenge  for  Systematic  Assessment  of  lonosphere/Thermosphere  Models  2:  Electron 
density,  Neutral  density,  NmF2,  and  hmF2  Using  Space  Based  Observations,  Space  Weather, 
doi:10.1029/2012SW000851 . 

Lei,  J.,  A.  G.  Burns,  J.  P.  Thayer,  W.  Wang,  M.  G.  Mlynczak,  L.  A.  Hunt,  X.  Dou,  and  E.  K.  Sutton 
(201 2),  Overcooling  in  the  upper  thermosphere  during  the  recovery  phase  of  the  2003  October 
storms,  J.  Geophys.  Res.,  117,  A03314,  doi:10. 1029/201 1JA016994. 

Huang,  C.  Y.,  S.  H.  Delay,  P.  A.  Roddy,  E.  K.  Sutton,  and  R.  Stoneback  (2012),  Longitudinal 
structures  in  the  equatorial  ionosphere  during  deep  solar  minimum,  J.  Atmos.  Sol. -Terr.  Phys., 
90-91 ,  doi:1 0.1 01 6/j.jastp.201 2.04.01 2. 

Wise,  J.  O.,  W.  J.  Burke,  and  E.  K.  Sutton  (2012),  Globally  averaged  exospheric  temperatures 
derived  from  CHAMP  and  GRACE  accelerometer  measurements,  J.  Geophys.  Res.,  117, 
A04312,  doi:1 0.1 029/201 1 JA01 71 08. 


A.5  Presentations 

Sutton,  E.  K.,  and  F.  Schmidt,  First  Principle  Modeling  of  the  Upper  Thermosphere/Lower  Exo¬ 
sphere  (invited),  CEDAR  Workshop,  Seattle,  WA,  21-25  June,  2015. 

Sutton,  E.  K.  (2015),  Satellite  Accelerometer  Observations  During  the  Storms  and  Substorms 
Without  Borders  Campaigns  (invited),  CEDAR  Workshop,  Seattle,  WA,  21-25  June,  2015. 

Thayer,  J.  P,  X.  Liu,  W.  Wang,  A.  Burns,  E.  K.  Sutton,  Composition  change  and  its  effect  on  ther¬ 
mosphere  mass  density  response  during  geomagnetic  activity,  European  Geophysical  Union 
Gen.  Assemb.,  Vienna,  Austria,  27  April  -2  May,  2014. 
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Koller,  J.,. . . ,  E.  K.  Sutton,  et  al.,  IMPACT:  Integrated  Modeling  of  Perturbations  in  Atmospheres 
for  Conjunction  Tracking,  AGU  Fall  Meeting,  abstract  #SM52C-08,  Dec.  2013. 

Huang,  C.  Y.;  Y.  Su,  E.  K.  Sutton,  D.  R.  Weimer,  R.  Davidson,  Energy  coupling  during  the  August 
2011  magnetic  storm,  AGU  Fall  Meeting,  abstract  #SM41A-221 5,  Dec.  2013. 

Weimer,  D.  R.,  T.  Edwards,  E.  K.  Sutton,  and  E.  Doornbos,  Storm-time  and  quiescent  variations 
in  the  thermosphere  exhibited  in  maps  of  exospheric  temperatures,  AGU  Fall  Meeting,  abstract 
#SA33C-06,  Dec.  2013. 

Liu,  X.,  J.  P.  Thayer,  W.  Wang,  A.  G.  Burns,  and  E.  K.  Sutton,  Composition  change  and  its  effect 
on  mass  density  response  during  geomagnetic  storm,  AGU  Fall  Meeting,  abstract  #SA33C-03, 
Dec.  2013. 

Shi,  Y.,  E.  Zesta,  H.  K.  Connor,  Y.  Su,  E.  K.  Sutton,  C.  Y.  Huang,  D.  Ober,  C.  Christodoulou, 
High-latitude  Thermosphere  Response  to  Solar  Wind  Dynamic  Pressure  Enhancements,  AGU 
Fall  Meeting,  abstract  #SA31  A- 1968,  Dec.  2013. 

Shi,  Y.,  E.  Zesta,  E.  K.  Sutton,  J.  O.  Wise,  Y.  Su,  S.  H.  Delay,  Effect  of  Solar  Wind  Dynamic  Pres¬ 
sure  Enhancements  on  Thermosphere  Density,  AGU  Fall  Meeting,  abstract  #SA33A-2172, 
Dec.  2012. 

Vickers,  H.,  M.  J.  Kosch,  E.  K.  Sutton,  Y.  Ogawa,  C.  La  Hoz,  Thermospheric  atomic  oxygen 
density  estimates  using  the  EISCAT  Svalbard  Radar,  AGU  Fall  Meeting,  abstract  #SA31A- 
2143,  Dec.  2012. 

Huang,  C.  Y.,  Y.  Su,  E.  K.  Sutton,  D.  R.  Weimer,  G.  D.  Earle,  A.  J.  Gerrard,  D.  M.  Ober,  G.  R. 
Wilson,  R.  Davidson,  Ion-neutral  Coupling  During  the  August  201 1  Magnetic  Storm.  AGU  Fall 
Meeting,  abstract  #SA31A-21 37,  Dec.  2012. 

Emery,  B.  A.,...,  E.  K.  Sutton,  et  al.,  Climatology  Assessment  of  lonosphere/Thermosphere 
Models  in  Low  Solar  Flux  Conditions  for  the  CCMC  CEDAR  Challenge,  AGU  Fall  Meeting, 
abstract  #SA23A-21 52,  Dec.  2012. 

Sutton,  E.  K.,  C.  S.  Lin,  and  S.  B.  Cable,  Comparison  of  TIE-GCM  SigmaO/N2  Column  Ratios 
with  Measurements  of  the  SSUSI  Instrument,  AGU  Fall  Meeting,  abstract  #SA23A-2151,  Dec. 
2012. 

Cable,  S.  B.,  E.  K.  Sutton,  T.  Matsuo,  and  C.  S.  Lin,  Model  bias  characteristics  for  data  assimi¬ 
lation  in  the  thermosphere,  AGU  Fall  Meeting,  abstract  #SA21A-1 869,  Dec.  2011. 

Lin,C.S.;  Knipp,D.;  Marcos, F.;  Cable,  S.  B.;  Sutton,  E.  K.,  Survey  of  Thermosphere  Joule  Heat¬ 
ing  Locations  and  Magnetic  Disturbance  Conditions,  AGU  Fall  Meeting,  abstract  #SA41A- 
1824,  Dec.  2011. 

Huang,  C.  Y.,  S.  H.  Delay,  P.  A.  Roddy,  E.  K.  Sutton,  R.  Stoneback,  and  R.  F.  Pfaff,  Ionosphere- 
thermosphere  coupling  during  deep  solar  minimum,  AGU  Fall  Meeting,  abstract  #SA21B-04, 
Dec.  2011. 

Shim,  J.,. . . ,  E.  K.  Sutton,  et  al.,  Systematic  Assessment  of  lonosphere/Thermosphere  Models 
during  the  2006  AGU  Storm,  AGU  Fall  Meeting,  abstract  #SA41  B-1 857,  Dec.  2011. 

Sutton,  E.  K.,  S.  B.  Cable,  F.  A.  Marcos,  and  C.  S.  Lin,  Improved  basis  functions  for  dynamic  cal¬ 
ibration  of  semi-empirical  thermospheric  models,  AGU  Fall  Meeting,  abstract  #SM31B-2096, 
Dec.  2011. 

Wise,  J.  O.,  E.  K.  Sutton,  F.  A.  Marcos,  and  W.  J.  Burke,  Modeled  Compostion  of  the  Thermo¬ 
sphere  at  Solar  Minimum  using  CHAMP  and  GRACE,  AGU  Fall  Meeting,  abstract  #SA1  IB- 
1593,  Dec.  2011. 
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A.6  Professional  Activities 


Journal  reviews:  GRL,  JGR,  Space  Weather,  ASR,  AIAA,  AAS 
Proposal  reviews:  NSF,  AFRL,  AFOSR,  SBIR/STTR 

Student  mentoring:  Lt.  Frederick  Schmidt,  AFIT,  masters  thesis  defense:  Feb.  2016;  Mr. 
Justin  Bock,  U.  Minnesota/AFRL  Space  Scholar,  2014;  Regular  interactions  with  NRC  postdocs 
at  AFRL 

A.7  Technology  assists,  transitions,  or  transfers 
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Appendix  B:  Time-dependent  thermospheric  composition 
for  N  components 

B.1  Introduction  and  Nomenclature 

In  this  section,  an  equation  describing  the  evolution  of  major  species  composition  in  a  log- 
pressure  coordinate  frame  is  derived  by  combining  the  species-dependent  continuity  and  diffu¬ 
sion  equations.  The  derivation  closely  follows  that  of  Dickinson  and  Ridley  [7J;  however,  addi¬ 
tional  terms  describing  time  dependence,  eddy  and  thermal  diffusion  are  included  to  reflect  the 
current  implementation  within  the  TIE-GCM.  We  also  deviate  slightly  from  their  treatment  to  high¬ 
light  several  equations  that  are  useful  in  tracking  species-dependent  as  well  as  mass-averaged 
transport.  The  following  definitions  are  used: 

Dy  mutual  diffusion  coefficient  of  ith  and  jth  components 
g  gravitational  acceleration 
Hi  scale  height  of  zth  component  [=  kT/(mig )] 

H  scale  height  of  mixture  [=  kT/(mg)\ 

He,  Ke  eddy  diffusion  coefficients 
k  Boltzmann  constant 

L  differential  matrix  operator  of  normalized  pressure  forces 
mi  molecular  mass  of  zth  component 
m  mean  molecular  mass  [=  (Y^=inimi) / n\ 
rii  number  density  of  zth  component 
n  total  number  density  [=  Y^iLi n*] 

Pi  partial  pressure  of  zth  component  [=  riikT ] 
p0  reference  pressure 
p  pressure 

Si  source  or  sink  for  number  density  of  zth  component 
s  vector  containing  the  first  (N - 1)  components  of  miSi/p 
T  temperature 

V  horizontal  component  of  the  momentum-weighted  mean  velocity 
w  vertical  component  of  the  momentum-weighted  mean  velocity  [=  Dz/Dt \ 

Wi  deviation  of  vertical  velocity  of  zth  component  from  mean  velocity 
w\  contribution  to  Wi  from  molecular  diffusion 
w"  contribution  to  from  eddy  diffusion 
w  vector  containing  the  first  (N - 1)  components  of  nirruwi 
w'  vector  containing  the  first  (N - 1)  components  of  rnrrnw- 
w"  vector  containing  the  first  (N - 1)  components  of  riimiW " 
z  vertical  spatial  coordinate 
z  vertical  log-pressure  coordinate  [=  In {po/p)\ 
ct  diffusion  matrix 

an  thermal  diffusion  coefficient  of  zth  component 
9  latitude 
A  longitude 
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Vi  volume  mixing  ratio  of  ith  component  [=  rii/n] 
p  mass  density  of  mixture  [=  Yh=i 
'ipi  relative  density  of  ith  component  [=  ramij  p\ 

\1>  vector  containing  the  first  (TV- 1)  components  of  ^ 
u  vertical  motion  relative  to  log-pressure  coordinates  [=  Dz/Dt ] 


B.2  Mass  Continuity 


Neglecting  horizontal  diffusion,  each  component  satisfies  the  following  continuity  equation: 


d_ 

dz 


(niiriiWi) 


d 

mtSi  -  —  ( niirii )  -  V  ■  (n;m?;V) 


d 


—  (n.imjW 

dz 


The  right-hand  side  of  (B-1)  can  be  written  in  terms  of  the  relative  densities: 


d_ 

dz 


(niiriiWi) 


rriiSi 


('tpip)  +  V  •  (-0ipV)  + 


(B-1) 


(B-2) 


We  wish  to  transform  Eq.  (B-2)  from  a  spatial  to  a  log-pressure  vertical  coordinate 
system  under  the  assumption  of  hydrostatic  equilibrium  using  the  following  relationship: 


dz  =  Hdz  (B-3) 

When  applying  this  transformation  to  partial  derivatives  with  respect  to  time  and  horizontal 
spatial  coordinates,  the  vertical  coordinate  being  held  constant  must  be  considered.  The  follow¬ 
ing  equations,  which  also  require  the  assumption  of  hydrostatic  equilibrium,  are  used  to  complete 
this  transformation  [[24L-  Eqs.  (3.6)  and  (3.17)]: 


i  (dr 
H  \dt/ 

(B-4) 

> 

T(v,.,v) 

(B-5) 

where  the  subscripts  £  and  z  refer  to  the  vertical  coordinate  being  held  constant  under  partial  dif¬ 
ferentiation.  Additionally,  the  relationship  between  the  spatial  and  log-pressure  vertical  velocities 
is  as  follows  [[24Lt  Eq.  (3.12)]: 


w  =  loH  + 


+  V  ■  V Zz 


(B-6) 


Making  the  appropriate  substitutions,  noting  that  the  equation  of  state  and  our  assumption  of 
hydrostatic  equilibrium  imply: 


pH  =  — e-2  (B-7) 

9 

and  dropping  the  subscript  ‘z  from  derivatives  taken  with  respect  to  time  and  horizontal  spatial 
coordinates,  Eq.[fE^2)  becomes: 
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The  definition  of  wt  implies: 


N 

^^riimiWi  =  0  (B-9) 

1=1 

Mass  sources  are  assumed  to  arise  solely  from  the  dissociation  of  one  molecule  into  others 
so  that: 


N 

^^miSi  =  0  (B-10) 

i= 1 

Relative  densities  ^  are  defined  so  that: 

N 

X>  =  1  (B-1 1) 

i=l 

By  combining  (B-8)  for  each  component  and  notin^I(i-$QB-10),Sfld  (B-1 1),  the  continuity 
equation  describing  the  total  fluid  in  log-pressure  coordinates  is  obtained: 

V- V  +  e2^-(e-2w)  =0  (B-1 2) 

dz 

Thus,  by  invoking  the  assumption  of  hydrostatic  equilibrium  and  adopting  pressure  coordi¬ 
nates,  the  mass  flow  of  the  fluid  appears  incompressible,  transforming  the  mass  continuity  equa¬ 
tion  from  a  prognostic  to  a  diagnostic  equation  (i.e.  no  time  derivatives  appear  in  the  equation). 

Using  Eq.  (1312),  the  divergence  terms  of  EQB-8)  can  be  simplified  in  favor  of  advection 
terms,  yielding  the  following  equation: 


d_ 

dz 


(' riiiriiWi ) 


Poe 

9 


+  V  •  +  u) 


djij 

dz 


(B-13) 


Now  let  w  be  the  (TV- 1)  vector  with  components  s  the  (TV-1)  vector  with  components 

rriiSi/p,  and  \l>  the  (TV -1)  vector  with  elements  V’i-  Then  the  first  (TV- 1)  equations  of  (b53)  can 
be  written  in  vector  form  as: 


d  Po  -z  f  d \Er  _  j  T  d 

— w  = -— e  —  +  V  ■  +  a;— ^  -  s 

dz  g  \  dt  dz 


(B-14) 


B.3  Molecular  and  Thermal  Diffusion 

With  the  assumption  that  the  atmosphere  is  in  a  state  of  hydrostatic  equilibrium,  i.e.  dp/dz  = 
-pg,  the  equation  of  motion  for  the  Tth  component  of  an  TV-component  mixture  in  the  presence 
of  molecular  and  thermal  diffusion  [[3]E  Eqs.  (18.2,5)  and  (18.3,13)]  can  be  written: 
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(B-15) 


N 


E 


riirij 

nD„ 


w. 


( i  J_  i  0171  dT\ 
dz  Hi  T  dz  ) 


The  pressure  force  exerted  on  molecules  of  the  ith  component,  expressed  by  the  right-hand 
side  of  (B^15),  forces  these  molecules  to  flow  through  the  rest  of  the  mixture  in  balance  with 
collisional  drags  given  by  the  left-hand  side. 

Noting  the  partial  pressure  pi  =  p^Hm/mt,  [(E15)  becomes: 


n 


N 

E 


A 


LmjDy  v nimiw i> 


A 


mjDy 


(riimiw'i) 


~d_ 

dz 


111  dm 

H  Hi  m  dz 


&Ti  dT\ 


A 


(B-16) 


Eqs.  (B-9)  and[(B}-1 1) — noting  that  the  former  applies  to  ticked  quantities  as  well — are  now 
used  to  eliminate  w' and  t/wfrom  the  first  ( N -1)  equations[oT](B-16),  giving  for  the  fth  component: 


where 


N- 1 

AjimjfijWj) 

3= 1 


'  d 

(l  . 

1 

1  dm 

an 

dTY 

dz 

U 

H, 

m  dz 

T 

dz  ) 

CXy 


1 


mj^Dnsr 

1 


+ 


yN- 1 


1 


1 


mkDfc  mj^Dgsr 


A 


rn-NDiN 


A, 


3  =  1 

3  F  i 


(B-17) 


(B-18) 


and  mN  refers  to  the  molecular  mass  of  the  Mh  species. 

Now  let  a  be  the  (iV-1)  x  (iV-1)  matrix  with  elements  and  L  the  diagonal  matrix  of 
differential  operators  with  elements: 


"  d 

f1  . 

1 

1  dm 

an 

dTY 

dz 

u 

Hi 

m  dz 

T 

dz  ) 

(B-19) 


The  solution  of  the  nonsingular  system  of  Eqs.  (B-17)  can  now  be  expressed  in  matrix  form: 


w '  =  a  1L\Ef  (B-20) 

Following  Dickinson  and  Ridley  [fz|,  a  nondimensional  form  of  the  diffusion  matrix  6*  can  be 
derived  using  a  nondimensional  parameter  A  related  to  the  mutual  diffusion  coefficient  through: 


m^D 

mjDij 


(B-21) 


where  D  is  a  characteristic  diffusion  coefficient.  It  is  assumed  that  D  varies  with  pressure  and 
temperature  in  the  following  way: 


D 


Do 


(B-22) 


where  D0  =  0.2  is  the  characteristic  diffusion  coefficient  at  S.T.P.  ,  700=273  K,  _p00=106  g/(cm  ■  s2). 
The  parameter  a  ^ defined  by  Eq.[(B}18)  is  nondimensionalized  by  the  substitution 
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o.y  =  [m^nD)  a.ij 


(B-23) 


where  the  nondimensional  parameter  ai3  is  then 


OL{j 


—  (piN  +  Yk^i  (<Pik  —  (pin)  i’k 

(  0y  0iN  )  Pp,  ‘ 


J  =  1 
J  +  i 


(B-24) 


Additionally,  Eqs.  (B-3)  and  (B-7)  are  again  used  to  transform  the  vertical  coordinate  of 
the  right-hand-side  of  t|cp(B-20)  into  log-pressure  levels,  resulting  in: 


w  =  r 


rp  \  0.25 

_1  (  T00  \  Pom 
mNg 


ck~1L\I/ 


where 


Ly  djj 


"  d 

(  rrii  1  dm 

an 

dTY 

dz 

\  m  m  dz 

T 

dz  J 

t  is  a  characteristic  diffusion  timescale  defined  by: 


r  = 


P^Hl 

Poo  Do 


and  Ho  is  a  characteristic  scale  height: 


(B-25) 


(B-26) 


(B-27) 


Ho  = 


kTi 


oo 


mNg 


(B-28) 


B.4  Eddy  Diffusion 

In  an  atmosphere  dominated  by  a  single  constituent,  as  is  the  case  with  molecular  nitrogen 
in  the  lower  thermosphere,  eddy  diffusion  establishes  a  flow  which  acts  to  smooth  gradients  in 
the  volume  mixing  ratio  of  the  minor  constituents,  uu  as  follows  [33,t4(|: 


,,  r>  1  ^ 

W:  =  -AK-  — 


Vi  dz 

In  terms  of  mass  flow  rates  and  mixing  ratios,  Eq.  |B39)  becomes: 

„  -  T>  (  d  1  drh\ 

riiiriiWi  =  - nmhE  —  +  —  — % 

\oz  m  oz  ) 

Transforming  to  log-pressure  coordinates  and  writing  in  vector  form,  Eq.  (B-30)  becomes: 


(B-29) 


(B-30) 


w»  =  _?2 

g  \  az  m  oz 


(B-31) 


where  KE  =  KE/H 2 
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B.5  Composition  Equation 


Setting  the  total  species-dependent  mass  flux  w  =  w'+  w"and  combining  Eqs.  (B-2^)I<ind 
(ail)  to  eliminate  wfrom  Ecf(fe-14)  yields  the  composition  equation: 


_9 

dz 


mN 


-a. 


e 


—z 


(s 


at 


-  V- 


(B-32) 


In  the  current  TIE-GCM  implementation,  the  subscripting  order  of  the  major  neutral  species  is 
as  follows:  i  =  {02, 0.  He},  with  N2  chosen  to  be  the  Mh  species  due  in  part  to  the  assumptions 
stated  in  Section  lB74j. 
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